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WAVE  PROPAGATION  IN  AN ELASTIC 
MEDIA  WITH  APPLICATIONS  TO  SEISMOLOGY 

by 

Walter  Silva 

Department  of  Geology  and  Geophysics 
University  of  California,  Berkeley,  California,  U.S.A. 

ABSTRACT 


A  formulation  is  presented  which  incorporates 
linear  anelastic  attenuation  into  plane  layer  models 
in  an  exact  manner.  Several  examples  of  body  wave 
propagation  in  absorbing  media  are  presented.  Surface 
time  histories  are  compared  between  predicted  acceler¬ 
ation  records  using  the  plane-layered  model  and  data 
recorded  by  a  vertical  array.  Spectral  ratios  between 
the  surface  and  bedrock,  computed  for  the  horizontal 
components,  show  fair  agreement  with  the  model  pre¬ 
dictions.  In  particular,  the  importance  of  attenuation 
in  predicting  ground  motion  in  soils  is  demonstrated. 

It  is  further  shown  that  converted  waves  are  of  minor 
importance  while  lateral  propagation  can  be  significant. 

The  formulation  is  extended  to  both  Love  and 
Rayleigh  wave  propagation.  Eigenvalue  and  surface 
displacement  calculations  for  a  high  loss  soil  structure 
indicate  that  Rayleigh  waves  are  more  strongly  affected 
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dictions.  In  particular,  the  importance  of  attenuation 
in  predicting  ground  motion  in  soils  is  demonstrated. 

It  is  further  shown  that  converted  waves  are  of  minor 
importance  while  lateral  propagation  can  be  significant. 

The  formulation  is  extended  to  both  Love  and 
Rayleigh  wave  propagation.  Eigenvalue  and  surface 
displacement  calculations  for  a  high  loss  soil  structure 
indicate  that  Rayleigh  waves  are  more  strongly  affected 


by  attenuation  than  are  Love  waves.  Inverse  calculations 
for  upper  mantle  and  crustal  structures,  with  both 
synthetic  and  real  data,  reveal  a  significant  dependence 
of  surface  wave  attenuation  upon  the  velocity  structure. 
Use  of  this  information  can  greatly  aid  in  velocity 
inversions  and  demonstrates  the  incompatibility  of  con¬ 
ventionally  extracted  attenuation  data  with  respect  to 
phase  data.  Further  indications  are  that,  for  the 
overdetermined  case,  the  same  layering  (number  and 
thickness)  may  not  be  optimum  for  both  the  velocity  and 
attenuation  inversion  parameters.  In  addition,  the  same 
layering  is  probably  not  optimum  for  both  Love  and 
Rayleigh  wave  inversions. 
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INTRODUCTION 


Recently  there  has  been  a  greater  appreciation  of 
the  effects  of  anelastic  attenuation  in  the  earth.  The 
recent  upsurge  in  interest  is  primarily  due  to  a  recog¬ 
nition  that  the  velocity  dispersion  accompanying  atten¬ 
uation  can  be  significant  (Liu  and  Archambeau,  1976). 

In  view  of  this  much  effort  has  been  directed  towards 
dispersion  corrections  to  free  oscillation  data  which 
is  then  inverted  to  obtain  corrected  earth  models 
(Hart  et  al. ,  1977).  The  results  have  been  very 
encouraging  in  that  body  wave  and  free  oscillation  earth 
models  have  largely  been  reconciled.  In  light  of  this 
realization  it  now  appears  that  attenuation  will  receive 
much  more  consideration.  It  is  the  effort  of  this 
thesis  that  this  consideration  be  in  terms  of  exact, 
rather  than  approximate,  theory. 

In  Chapter  1  the  Haskell-Thompson  propagation 
matrix  technique  is  extended  to  include  anelastic 
attenuation  in  an  exact  manner.  In  order  to  demon¬ 
strate  the  effect  of  attenuation  (as  compared  to  purely 
elastic  propagation)  on  body  waves,  several  examples 
are  presented.  In  particular,  transfer  functions  are 
calculated  for  typical  soil  and  crustal  structures. 

In  addition,  an  example  is  shown  demonstrating  the 
effect  of  an  attenuating  boundary  layer  above  the 
core-mantle  boundary  on  reflected  pulses. 


V 

As  a  means  of  estimating  the  suitability  of  the 
plane-layer  model  in  predicting  ground  motion  compar¬ 
isons  are  made  between  the  predicted  surface  motion  and 
data  recorded  at  a  vertical  array.  The  array  is  located 
in  a  soil  section  and  demonstrates  the  effect  of  low 
velocity  surficial  material  on  wave  motion. 

In  Chapter  2  the  formulation  is  applied  to  both 
Love  and  Rayleigh  wave  propagation.  In  this  section  a 
soil  and  an  upper  mantle  model  are  considered  which 
demonstrate  the  effects  of  attenuation  on  surface  wave 
propagation.  In  addition,  an  inversion  scheme  is 
presented  by  which  depth  dependent  velocity  and 
attenuation  may  be  estimated  simultaneously  from  sur¬ 
face  wave  phase  and  amplitude  data.  The  inversion 
scheme  is  demonstrated  using  synthetic  data  and  is  then 
applied  to  real  data  suitable  for  upper  mantle  structures. 
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CHAPTER  1 
BODY  WAVES 


I.  FORMULATION 

In  both  the  body  wave  and  surface  wave  analysis 
the  matrizant  technique  (see  Haskell,  1953,  for  first 
seismological  application)  is  employed  as  the  compu¬ 
tational  algorithm.  In  this  approach  the  medium 
properties  are  assumed  piecewise  constant.  The  equations 
of  motion  are  integrated  analytically  and  the  solutions 
propagated  by  matching  boundary  conditions  at  the  layer 
interfaces.  The  main  disadvantage  of  this  method  over 
direct  numerical  integration  of  the  equations  of  motion 
(Gilbert  and  Backus,  1966)  is,  of  course,  that  it  is 
possible  only  for  plane  geometry.  To  overcome  this 
shortcoming,  various  earth  stretching  transformations 
have  been  developed  (Schwab  and  Knopoff,  1972)  which  are 
used  in  this  treatise  where  appropriate.  The  following 
paper  extends  the  matrizant  method  to  include  anelastic 
attenuation  in  an  exact  manner  and  forms  the  basis  for 
this  thesis. 
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BODY  WAVES  IN  A  LAYERED  ANHLASTIC'  SOLID 
Bn  W  Sii  v  a 

AHSI  K  VI  I 

A  formulation  extending  the  ttaskell-Thompson  matrix  method  to  include  the 
effects  of  anelastic  attenuation  is  presented.  The  formulation  is  exact  in  that  no 
low-loss  approximations  are  made.  Consideration  is  given  to  nonparallel 
propagation  and  attenuation  directions  with  corresponding  velocity  anisotropy . 
Examples  are  presented  for  models  representing  soils,  the  crust,  and  the  core- 
mantle  boundary . 


Introih  chon 

With  the  increase  in  the  number  of  stations  and  the  higher  degree  of  standardiza 
recent  years,  more  use  is  being  made  of  seismic  amplitude  data.  This  has  contribute 
increased  regionalization  of  structure  down  to  the  core-nutntle  boundary.  In  ore  .> 
accurately  represent  this  fine  si ;  ucturc  in  applying  corrections  or  to  resolve  it  in  inverting 
data,  more  use  is  being  made  of  the  higher  frequencies  where  the  attenuation  effects  are 
most  significant.  It  is  therefore  becoming  increasingly  important  to  consider  nongeometri- 
cal  attenuation  exactly.  Past  approximations  in  dealing  with  loss  (Knopoff,  1964)  must  be 
replaced  with  exact  formulations  (Lockett,  1 962; Cooper,  1967;  Borcherdt.  1971 ;  Buchen. 
1971). 

In  order  to  consider  the  effects  of  a  vertical  variation  in  attenuation  as  well  as  velocity 
and  density  on  body  waves,  an  extension  of  the  Haskell-Thompson  (Haskell.  1953)  matrix 
formulation  using  an  exact  theory  is  presented.  In  particular,  the  restricted  problem  of 
anelastic  layers  on  an  elastic  half-space  is  considered,  but  the  formulation  can  easily  be 
extended  to  include  an  attenuating  half-space  Previous  consideration  of  the  problem 
(Kanai.  1 950)  dealt  with  normally  incident  homogeneous  waves  with  \  iscoelasticity  of  the 
Voigt  type.  The  present  treatment  considers  incident  P  or  SI  waves  at  arbitrary  angles 
and  a  general  constitutive  relation. 


Formi  I  ATION 

The  most  general  form  of  a  linear  constitutive  relation  is  Boltzman’s  superposition 
principle  (Gurtin  and  Sternberg.  1962)  which,  w  ritten  in  terms  or  the  tensorial  relaxation 
function  r(f )  is 

,  0) 

=  r,,t,(t  |*,/lu(M 

where  P,y(( )  and  i:,y(t )  are  the  lime-dependent  stress  and  strain  tensors  and  the  symbol  * 
denotes  the  Stieltjcs  convolution. 

Assuming  the  medium  to  be  isotropic  and  homogeneous,  equation  (1 )  may  be  broken 
up  into  bulk  and  shear  components  and  written  as 

P,j[l  I  =  2/t(f  )  *  I  I  *  / 

P„(M  =  3MM*./«:H(M  (2) 

where  //(f)  and  Mf)  are  the  relaxation  functions  in  shear  and  bulk.  Assuming  that  the 

is. w 


Utmp  — 


/ 


1 540 


W  Sll  VA 


particle  displacements  »,  are  infinitesimal,  the  strain  can  be  written 
and,  neglecting  body  forces,  the  linear  momentum  equation  is 

(3) 

where  />  is  the  medium  density.  Substituting  equation  (2)  into  equation  (3)  yields  the 
equation  of  motion. 

[MO  +  j/i(tl]*[V(V  du)]  [ , iff ).[Vx(Vxdu )]=,,«.  (4 » 

Since  the  convolutions  make  the  time-domain  representation  quite  intractable,  it  is 
customary  to  take  the  Fourier  transform  of  equation  (4).  Restated  in  terms  of  transformed 
variables,  equation  (4)  becomes 

[k  +  t,«]V(V  ■  u)  -  [//]V  x  (V  x  u)  =-  -  pw2u  (5 ) 

where 

K  =  i(nj’  ,  k(i)c  "" dt ,  ft  =  i<uf‘  ,  /i(t)e  ""dt 
u  =  J  ‘  ,  u 

At  this  point  it  is  convenient  to  introduce  the  transformed  P  and  S  displacement 
potentials  in  terms  of  Helmholtz's  (elation 


u-V$  +  V  x  tji.  Viji  =  0. 


(6) 


Substituting  equation  (6)  into  equation  (5)  results  in  the  familiar  Helmholtz  equations 
for  the  P  and  S  potentials  <f>  and 

[V2+Kp2](£  =  0.  [V2+*\2]^  =  0  (7) 

where 


KP2 

Ks2 


=  (!> 


2 


=  t/t 


2 


P _ 

K(t0)  +  *fi(<0) 
P 

fiU») 


(8) 


Note  that  the  terms  i2  and  ft2  are  in  general  frequency-dependent  in  both  real  and 
imaginary  parts. 


Mfdiijm  Paramftfrization 

Let  us  now  consider,  for  demonstration  purposes,  the  case  of  S  waves.  A  general 
solution  for  tjr  in  equation  (7|  is 

^  =  ^(m)expt-iKsX)  (9) 

where  Ks  is  a  complex  vector  with  the  real  and  imaginary  parts  having  diflercni  directions 
in  general. 

Ks  —  Ps  —  /,  Av 

Ks2  =  KyKs  =  |Ps|2-|zls|2-/2PsAs 
P,  •  A.v  =  |PV|  |AS|  cos  (•/,). 


(10) 

(11) 

(12) 
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Pv  is  the  propagation  vector  such  that  <•>  |l’vj  is  the  phase  velocity  and  As  the  attenuation 
vector  such  (hat  exp  (  As  \l  represents  the  spatial  decay  of  the  potential.  The  non/ero 
gives  rise  to  the  inhomogeneous  waves  (Borcherdt  1971.  Buchen  1971.  Cooper  1967. 
Lockett  1962)  whose  amplitude  varies  (monotonically )  along  a  wave  front.  It  becomes 
necessary  now  to  specify  the  three  parameters  | l*s] . |AS|.  and  in  terms  of  material 
properties  and  medium  geometry. 

Writing  the  transformed  shear  modulus  //(o)  in  equation  (X)  in  terms  of  a  real  part, 
and  an  imaginary  part.  /i,(r>).  the  quality  factor  Qs  for  shear  waves  is  defined  as 


Q  s 


«,(<•»)  I  At. 

I  277  I. 


(12) 


where  is  the  peak  energy  density  stored  and  Al:  is  the  energy  lost,  both  per  cycle 
(Borcherdt.  |97|.  1977).  K  s  ’ may  be  written  in  the  following  form 


fs\- 


I  *  s  I  *  Qs 


(14) 


where  i\  is  the  homogeneous  wave  velocity  of  the  medium.  I  sing  (14)  ti>  invert  |l  I  land 
( 1 2)  we  arrive  at  convenient  expressions  for  T\j  and  As' 


l\ 

1  '  x  1 

1 

(  1  I  .  1  t  (>s  COS  -  (  J) 

-  Qs 

(15) 

Ay’ 

1  •»* 

1 

1  ‘  x  1 

(  I  1  \  1  1  Qs  '  COS  •  ( )  ) 

*  Qs  • 

(16) 

with  similar  expressions  for  /’  waves  using  the  /’  parameters  In  the  low -loss 
approximation  for  homogeneous  waves  (  .  0.  (f.  *•  I  )  equation  1 16)  reduces  to  the  well- 

know  n  expression 


When  dealing  with  highly  dissipative  materials  the  vectorial  nature  of  A  must  be 
considered.  The  problem  is  that  for  a  given  incident  wave  (direction  of  both  I’  and  A 
specified )  onto  a  plane  boundary  between  two  v iscoelastic  media,  the  direction  of  both  /’ 
and  1  must  be  determined  for  the  /’  and  \1  reflected  and  transmitted  waves.  These 
directions  can  be  uniquely  determined  by  applying  the  usual  boundary  conditions  at  a 
welded  interface  (or  free  surface  for  half-space  problems).  This  results  in  an  extended  form 
of  Snells  law  in  that  I ,  as  well  as  l\  must  be  continuous  (Borcherdt  1977.  Lockett  1962 1 
In  the  restricted  case  considered  in  this  paper  where  the  incident  medium  is  elastic.  -I,  is 
zero  everywhere.  This  enables  elastic  layers  to  be  inlerbedded  with  absorbing  media. 


Lxiivsion  ni  mi  Hasm  1 1  - Thompson  Loknii  i  a  i  ion  i  ok  I.vMRll)  Midi  a 

The  following  development  follows  closely  that  of  Haskell  (1953,  1962).  However, 
displacement  potentials  are  used  here  Referring  to  Figure  I  for  coordinate  reference,  we 
can  write  the  solutions  to  equation  (7)  in  the  usual  form 


</>  |  .-IpCxplifc,../ )  f  Bp  exp  (  iA  ,../)]  exp  |  -  iAr>.\  ) 

iji  ^|r,  |  lsexp(iA\./t  r  Bvexp(  iA\./l|exp(  iA\,V) 


(17) 
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where  A,  B  are  complex  and  in  general  frequency-dependent  amplitudes.  From  equation 
(10) 

Kz  =  P.  —  iaz\  Kx  =  Px  —  iAx. 

In  the  simple  case  we  are  considering  (incident  elastic  wave)  we  have  Ax  =  0,  and  from  the 
boundary  conditions  .4,  must  be  zero  everywhere.  Thus  Kx  remains  a  real  quantity.  If  the 
incident  medium  were  anelastic,  the  incident  attenuation  direction  would  have  to  be 
specified  along  with  the  propagation  direction,  and  then  A.  in  each  layer  would  adjust 
itself  to  be  consistent  with  a  continuous  Ax  and  a  specified  Kp  s  for  that  layer  (see  equation 
14).  Choosing  AP  and  ,4S  as  upgoing  potentials  (negative  ;  direction)  we  can  write  for  each 
layer 

K.  =  principal  value  (K1  -A.',2)1  2.  (18) 


Fto  I  The  problem  is  uniquely  specified  given  I-,,  Vt.  Qr.  Qs,  p,  and  Z  for  each  layer  and  given  Ar 5  (Incidenl 
P-  or  S-wave  potential  amplitudel  in  the  elastic  half-space  For  an  anelastic  half-space,  the  direciion  of  the 
incident-wave  attenuation  vector  must  also  be  specified. 


Using  equations  (2),  (3),  (6),  (17),  the  displacements  and  stresses  for  layer  mean  be  put  in 
the  following  matrix  form 


u  1 

i  —  ‘KPxCP 

Ks,Ss 

KPxSP 

-«s2Cs 

M' 

1  ~Kp;Sp 

iKPjCP 

Ks*Ss 

P.. 

i2pK  rJtp:Sr 

~(inCs 

2fiKPJCPxCP 

—  iuilSs 

P"L 

i2fiKS!KPxSs 

ifiSiSp 

2  nKSlKSxCs 

4P+  BP 

As  + 

Ap—  Bp 

As  ~ 

Q  =  Cos  (Ks.Zm) 


where 


Cp=Cos<Kp,ZJ 
Sp  =  Sm(Kp,Zm) 
S1  =  K2X-K2S; 


Ss  =  S,n(KS!ZJ. 
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This  result  (which  is  equivalent  to  equation  (3  20)  of  Grant  and  West,  1965)  can  be 
conveniently  written  as 


■V.-OJ/.K'..  (19) 

Thus  we  see  that  (layer  thickness)  is  the  phase  factor  w  hich  propagates  the  potentials 
across  the  mth  layer  and  that  l)m  may  be  thought  of  as  a  form  of  propagator  matrix  w  ith 
C'm  the  coefficient  matrix.  With  this  in  mind  and  w  ith  the  idea  of  eliminating  C'm  we  can 
write  (Haskell.  1953) 

•V-  -  Dm  1  (0TVm  (20) 

Then  applying  the  usual  boundary  conditions 

■Vm  =  ( Dj/m)Dm  ‘(0 ))(Dm  ,(/„  ,)/)„',( oiia;  , 

~  i . . .  (i,A'n.  (21) 


and  for  n  -  1  layers  where  layer  n  is  an  elastic  half-space  and  interface  0  is  a  free  surface 

1  (0 )o„  i <Jn  :■  a,Ao 

=  J  X„ 

with  the  following  matrix  elements. 

0  0 
-i2KFx  -1  //i 
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The  elements  of 
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o 
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(22) 


(23) 


(24) 


are  given  by 


<»,,  =llCs-2A^Cp 

</,2  =  —  *Ap,f 2AS.5X  +  (d  A p. ) 5 r] 

u,  5  ~  —  ft  1  [A s.Ss  +  (Kpx  Kp.)Sp] 


‘*14  -  tAp./t  [kp  f  s  I 

11 2 1  =  ,Ap1[2Ap.Sp  +  |fl  As.  )Ss] 

«22=nCp-2A^Cs. 

<j  2 ,  =  <r ,  4 


r/24  —  /i  ‘[Ap..5p  +  (Kpx  Ks.  )Ss] 
~  /1[4A  pjA  p.Sp  +  (ft2  ASj  )Ss] 

« _i  j  =  —  /2/tA  p,fl[  Gp  —  Cs] 
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l,3i  —  41 1  1 
"34  =  "21 


"4 1  —  ui2 

"42  =  /i[4K^SsS,  +  <Q2  KP;  IS,] 

"43  =11 1  j 
"44  =  "2  2 

C„  therefore  becomes  the  input  matrix  and  choosing  AP  s  in  the  upgoing  ( -  c)direction 
and  considering  incident  P,  we  can  invert  equation  (22)  to  give  the  surface  displacements 
t/0  and  w'o  in  terms  of  the  incident  potential  [„Ap),Kx.  and  the  layering. 

"o  =  ~  “[*^22  4“  J  f!  R 

vv0  =  ~[^2 I  +  4  I  ]n^p  R 


R  —[J 2  1  +  J4|][*f  12  +  ^32]  IV 22  ^  ^42] IV  II  + 


(26) 


APPLICATIONS 

In  order  to  illustrate  the  effects  of  attenuation,  three  models  which  represent  soils,  the 
crust,  and  the  core-mantle  boundary  are  considered.  The  structures  arc  listed  in  Table  I. 
With  the  exception  of  the  low-velocity  layer  of  the  upper  mantle,  these  appear  to  be  the 
three  regions  where  nongeometrical  attenuation  is  most  pronounced  and  therefore  may 
have  some  effect  on  observational  interpretation.  Also,  knowledge  of  the  Q  structure  of 
Ihese  regions  will  be  valuable  in  interpreting  materials  and  structure  mechanisms  when  an 
acceptable  theory  is  found  relating  state  variables,  material  properties,  and  energy 
absorption. 

In  applying  this  formulation  in  calculating  reflection  and  transmission  coefficients, 
transfer  ratios,  synthetic  seismograms,  etc.,  some  estimation  must  be  made  of  the  medium 
parameters.  This  usually  means  a  frequency-independent  loss  and  velocity  which  can  be 
shown  to  violate  causality  (Futtcrman,  1962).  However,  since  the  frequency-dependence 
can  be  made  weak  over  a  finite  frequency  band,  assuming  a  frequency-independent  loss 
and  phase  velocity  over  the  space-time  dimensions  considered  here  should  not  be  critical. 

(u)  Soils.  The  effects  of  attenuation  can  be  rather  drastic  in  a  highly  dissipative  material 
such  as  loosely  compacted  soils.  The  structure  chosen  (Table  1 )  is  for  the  Richmond  Field 
Station  of  the  University  of  California.  Berkeley  and  consists  of  mud  deposited  in  San 
Francisco  Bay.  Borehole  measurements  of  velocity  and  sample  measurements  of  both 
velocity  and  density  were  available  for  this  site.  The  Q  structure  represents  a  best  guess  for 
illustrative  purposes  (structuredata  from  1 .  V.  McEvilly,  oral  comm.).  Figure  2  shows  the 
vertical  and  horizontal  displacement  spectra  for  normally  incident  P  and  S  waves, 
respectively.  All  input  potentials  were  normalized  to  unity  total  displacement  for  incident 
1 )  or  SF(„4s  =  „ks  1 ).  The  solid  line  is  for  an  elastic  stack  while  the  broken 
line  includes  the  effect  of  loss.  The  vertical  motion  is  somewhat  unstructured  because  the 
compressional  wavelengths  are  greater  than  any  of  the  layer  thicknesses.  The  loss  behaves 
as  we  might  expect  for  purely  homogeneous  waves,  mirroring  the  clastic  behavior  at  a 
lower  amplitude  and  becoming  asymptotic  to  it  toward  low  frequencies.  Considering  the 
shear  spectra  (Figure  2B)  wc  begin  to  note  some  interesting  effects.  First,  the  elastic 
spectrum  shows  the  characteristic  peaks  (shear  wavelengths  <  layer  thickness)  which  are 
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TABLE  I 

Physicai  Parami  ors  k>r  Ri  RRisisiATivi  Modi  is 
CONSIDERID  IN  NUMERICAL  CAU'll  ATIONS 
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■/ 
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8.08 

0 

9.90 

2500 

, 

resonances  associated  with  the  total  S-wave  travel  time  (Bakun.  1971,  Haskell.  1960).  The 
total  S-  wave  travel  time  of  the  stack  is  T  =  0.098  sec  and  maxima  and  minima  arc  expected 
at 


n  m 

/ma»  =  47- .  n  =  1 . 3. 5 — ;  /mi„  =  27..»«= 1.2.3.... 

=  2.6.  7.9. 12.7...  =5.1. 10.2,20.3. .. . 

The  peaks  and  troughs  are  not  exact  because  the  total  stack  travel-time  effect  is  modulated 
by  the  layering.  In  the  loss  spectra  we  see  that  there  is  little  information  content  at 
frequencies  greater  than  about  14  Hz.  The  effect  of  attenuation  is  more  drastic  for  shear 
waves  due  to  the  lower  Qs  and  the  longer  travel  times.  Also,  it  is  important  to  note  the 
slight  shifting  of  the  peaks  in  the  case  of  loss.  The  velocities  arc  the  same  in  the  clastic  and 
attenuating  layers  and  the  shifting  is  due  to  the  change  in  modulation  as  the  loss  affects  the 
acoustic  impedance. 

In  Figure  3  are  shown  the  crustal  transfer  function  ratio  (»v„  n„ ).  the  vertical  spectra, 
and  the  horizontal  spectra  for  an  incident  comprcssional  wave  at  i  =  10  for  the  same  soil 
structure.  It  is  interesting  to  note  the  considerable  change  in  the  ratio  for  the  loss.  Any 
inversion  scheme  not  accounting  for  the  loss  would  yield  a  different  structure.  Again  the 
shear  spectrum  is  the  controlling  mechanism  but  in  this  case  the  large  discrepancy  between 
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the  elastic  and  loss  is  largely  due  to  the  velocity  anisotropy  induced  by  the  inhomogeneous 
waves  (see  equation  15). 

(6)  The  crust.  The  crustal  model  (Table  1 ).  excepting  the  Q  structure,  was  taken  from 
Bakun's  best-fitting  Berkeley  crustal  model  (Bakun,  1970).  The  Q  structure  represents  a 


Fig.  2  Spectra  of  normalized  surface  displacements  for  a  high-loss  soil  structure  (Table  1).  Solid  lines  are 
elastic  layers,  broken  lines  include  loss  (At  Vertical  displacement  for  incident  P  wave,  (B)  Horizontal 
displacement  for  incident  S  wave;  both  at  normal  incidence. 

best  guess  for  Qs  by  the  author  based  on  some  near-Berkeley  crustal  studies  (Kurita,  1975, 
O’Neill  and  Healy,  1973)  and  the  relation 

e„=2<?s<WJ  (27) 

The  transfer  ratio  for  the  crustal  model  for  an  incident  compressional  wave  at  i  =  25°  is 
shown  in  Figure  4  along  with  the  vertical  and  horizontal  spectra.  The  transfer  ratio  for  the 
elastic  and  loss  agree  well  out  to  about  3  Hz  which  is  high  enough  to  resolve  the  structure. 

The  spectra  of  the  vertical  and  horizontal  surface  displacements  for  the  Berkeley  crust 
have  been  synthesized  and  are  shown  in  Figure  5,  where  (A)  and  (B)are  the  vertical  elastic 


Kill  4.  Spectra  of  normalized  |A|  crustal  transfer  function  (h0/u0).  IB)  vertical  surface  displacement.  (Cl 
horizontal  surface  displacement  for  incident  P  wave  at  25  for  the  Berkeley  crust  (Table  I ).  Solid  lines  are  elastic 
layers,  broken  lines  include  loss 
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and  loss  w  hile  (C|  and  (D|  represent  the  horizontal.  The  effect  of  the  loss  shows  a  general 
smoothing  of  the  record,  a  decrease  of  the  higher-order  reflections,  and  a  reduction  in 
amplitude  (vertical  peak  to  peak  by  0.7,  horizontal  by  0.8). 

(<•)  Core-mantle  boundary.  To  demonstrate  the  effects  of  loss  near  the  core-mantle 
boundary  on  PcPand  PeS.  the  reflection  spectrum  for  the  core-mantle  boundary  structure 
of  Table  I  was  synthesized.  The  velocities  are  from  Bolt  (1972)  and  the  boundary  layer 
density  was  derived  assuming  the  region  to  be  a  thermal  boundary  layer  and  to  consist 
entirely  of  mantle  material  (Glyn  Jones,  personal  comm.).  The  QP  is  from  Kuster  (1972) 
while  Qs  is  derived  from  equation  (27)  (zero  loss  in  bulk).  An  incident  P  wave  is  considered 
with  i  =  25  and  the  synthesis  is  performed  for  a  point  70km  from  the  boundary  layer. 
Figures  6  and  7 show  the  synthesized  potential  coefficients  („Bp  l,Ap.„BSnAp)  for  the  same 
explosion  source  as  the  crustal  seismograms.  The  first  small  pulse  in  Figures  6  and  7 
represents  reflected  P  and  SF  waves,  respectively,  from  the  abrupt  transition  between  the 
lower  mantle  and  the  boundary  layer.  A  more  realistic  gradient  would  largely  eliminate 
this  reflection.  The  figures  then  represent  PeP  and  PeS  sources  to  be  convolved  with 
suitable  transfer  functions  and  show  the  effect  of  attenuation  on  the  reflected  amplitudes 
(PeP  zero-to-peak  reduction  0.8,  PeS  0.6).  The  effect  on  the  wave  forms  seems  to  be  small 
at  this  angle  of  incidence  for  such  low  Q  values  and  indicates  that  a  considerable  amouni  of 
attenuation  is  possible  in  the  boundary  layer  and  still  be  unobservable. 


Apim  vdix 

To  consider  a  fluid  layer  (/<  =0)  the  matrix  D, (equation  19)  must  be  modified  due 
to  the  overspecified  boundary  conditions  at  a  solid-fluid  interface.  Using  a  development 
similar  to  Tcng  (1967)  )  for  a  fluid  layer  becomes 

0  0  0  1 

-  kP:SF  0  ‘kp.Cp  0 

0  10  0 
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Fl(.  5.  A  ANI)  B 

Fici.  5  Synthetic  seismograms  for  the  Berkeley  crust  |A|  and  (Blare  vertical  motion  (positive down) for  the 
elastic  and  loss  cases,  respectively,  which  were  synthesized  from  the  spectra  in  Figure  4B.  (C)  and  (D)  are 
horizontal  motion  for  the  elastic  and  loss  cases,  respectively,  which  were  synthesized  from  the  spectra  in  Figure 
4C  All  were  convolved  with  an  explosion  source  function  appropriate  for  BOXCAR  (Helmberger  and 
Harkrider.  1972).  a  Bcnioff  short -period  instrument,  and  a  low-pass  filter  with  a  corner  frequency  at  5  Hz. 
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Fl(i  h  Synthesized  reflected  P  potential  coefficient  for  incident  P  reave  on  core-mantle  boundarv  structure  of 
Table  I  The  frequency  interval  was  001  Hz  with  M2  sample  points  All  spectra  were  filtered  with  a  low -pa  s 
causal  filter  with  corner  frequency  at  2  0  Hz  I  A)  clastic,  |B|  loss 
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H<*  7  S>nthcM/ed  reflected  S  potential  coefficient  for  incident  /’  wave  on  core-mantle  boundary  structuie  of 
Table  I  The  frequency  internal  was  0 01  H/  with  51^  sample  points  All  spectra  were  filtered  with  a  low-pass 
causal  filter  w ith corner  frequency  at  2  OH/  t A l  clastic.  (B| loss 
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II.  APPLICATION 

In  this  section  the  body  wave  formulation  will  be 
applied  to  an  earthquake  (Briones  Hills:  8  Jan.  1977 » 

M  =  4.3,  37°  54.31 'N  122°  10.97 'W,  depth  9*5  km, 

A  =  13  km)  recorded  at  a  vertical  array.  The  array 
consists  of  three  3-component  accelerometers  positioned 
at  depths  of  120  feet  (bedrock),  40  feet  (bay  mud),  and 
the  surface  (soil)  located  at  the  Richmond  Field 
Stations.  The  array  is  near  the  borehole  referred  to 
in  Secion  I  so  the  soil  parameters  listed  in  Table  I 
(Section  I)  can  be  used  to  calculate  transfer  functions. 
The  general  problem  of  calculating  the  effect  of  local 
geology  on  earthquake  ground  motion  is  critical  in  the 
design  of  structures.  Of  particular  concern  is  the 
amplification  factor  associated  with  low  velocity 
surficial  material.  In  order  to  test  the  suitability  of 
the  body  wave  formulation  in  predicting  the  response 
of  the  mud-soil  structure,  a  comparison  is  made  between 
the  observed  surface  acceleration  records  and  the  bedrock 
accelerograms  continued  to  the  surface.  Different 
transfer  functions  (eq.  26,  Section  I),  applied  in  the 
frequency  domain,  encompassing  different  angles  of 
incidence  will  be  considered.  Also  both  elastic  and 
anelastic  propagation  will  be  compared  to  the  observed 
surface  time  histories.  The  entire  accelerogram  is 
used  in  computing  the  spectra  and  no  time  window  has 
been  applied.  Twenty  seconds  of  record  is  transformed 
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with  2048  points.  The  analysis  will  be  restricted  to 
horizontal  components  as  the  bedrock  vertical  instru¬ 
ment  had  a  telemetry  malfunction.  Also,  the  actual 
orientation  of  the  horizontal  components  have  not  yet 
been  determined  and  will  therefore  be  referred  to 
simply  as  HI  and  H2. 

In  Figure  1  are  shown  the  bottom  and  surface  HI 
(horizontal)  component  acceleration  records  plotted 
to  the  same  scale.  The  onset  of  the  P-wave  is  clearly 
visible  while  the  S  arrival  is  somewhat  emergent.  Of 
particular  interest  is  the  long  period  component  of 
the  bottom  S  arrival  compared  to  the  corresponding 
part  of  the  surface  record.  The  surface  trace  appears 
to  have  more  high  frequency  content.  This  observation 
is  substantiated  by  the  spectra  shown  at  the  end  of  the 
section.  Compare  the  H1B  and  HIT  spectrums.  The  peak 
in  the  HIT  between  1  -  2  Hz  is  missing  in  the  bottom 
spectra.  Also  the  bottom  record  tapers  off  after  the 
S-wave  while  the  surface  record  continues  with  large 
amplitudes  and  points  out  the  importance  of  high 
frequency  surface  waves  in  the  duration  of  motion 
for  near  sources.  The  surface  wave  traces,  either 
Love  or  Rayleigh  waves,  represents  a  laterally  prop¬ 
agating  disturbance  and  is  not  included  in  the  theo¬ 
retical  model.  A  full  treatment  of  the  problem  requires 
the  inclusion  of  a  non-plane  wave  source  and  is  not 
considered  here. 
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In  order  to  simulate  the  layering  effect,  the 
bottom  trace  (Figure  1,  a)  was  continued  to  the  surface 
using  the  previously  mentioned  transfer  function. 

Figure  2  shows  the  computed  surface  motion  assuming 
Figure  la  as  a  normally  incident  S-wave .  Figure  2a, 
considers  the  soil  elastic  while  Figure  2b  includes  the 
loss  (note  different  scales).  In  comparing  the  observed 
surface  motion  with  the  computed  (including  loss) 
several  interesting  features  are  observed.  First, 
the  amplitudes  are  fairly  close  (not  considering  the 
P-wave  as  it  was  continued  as  an  S-wave)  with  the 
calculated  somewhat  large.  This  indicates  that  the 
Q  structure  is  approximately  correct.  Also,  the  longer 
period  part  of  the  S-wave  is  present  on  the  continued 
trace  as  we  might  expect  since  its  wavelength  is  far 
greater  than  the  thickness  of  the  soil  section.  In 
addition,  the  large  amplitudes  following  the  S-wave 
are  not  present  on  the  computed  seismograms.  This 
indicates  that  this  part  of  the  motion  is  propagating 
in  a  horizontal  direction.  Figure  3  shows  the  same 
input  (assumed  SV)  continued  with  an  incidence  angle 
of  10°.  This  trace  demonstrated  that  converted  waves 
are  unimportant  in  this  case  as  it  is  virtually  the 
same  as  Figure  2. 

Considering  now  Figure  4,  we  have  the  same  analysis 
for  the  other  (H2)  component.  Again  (a)  and  (b)  arc 
t.he  bedrock  and  surface  records  respectively.  Again 
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the  longer  periods  seem  more  dominant  in  the  bedrock 
component.  Also,  the  surface  waves  are  much  more 
dominant  in  this  component,  appear  to  arrive  later, 
and  may  be  associated  with  Rayleigh  waves.  Comparing 
the  surface  record  (b)  with  the  computed  in  Figure  5 
we  observe  the  amplitudes  of  the  attenuation  record 
(b)  are  fairly  close.  The  longer  periods  are  continued 
as  dominant  while  the  large  amplitude  and  long  period 
wave  train  is  absent.  Figure  6  shows  the  continued 
record  for  a  transfer  function  with  an  incident  SV 
wave  at  10°  and  is  virtually  the  same  as  0°  in  Figure  5- 
In  order  to  compare  frequency  domain  data  with 
the  model,  spectral  ratios  (surface  to  bedrock)  were 
calculated  for  both  components  and  are  shown  in  Figure 
7.  The  ratios  shown  have  been  smoothed  with  a  20 
point  averaging  filter.  Because  of  this  smoothing, 
coupled  with  the  spectral  contamination  of  a  boxcar 
window,  the  absolute  magnitudes  cannot  be  directly 
compared  with  the  theoretical  transfer  function. 

However,  the  observed  ratios  (solid  line  in  Figure  7), 
except  for  the  peak  near  5  Hz,  agree  fairly  well  with 
the  theoretical  transfer  function  (broken  line  in 
Figure  7;  also  shown  in  Figure  2b,  Section  I  (an 
unfortunate  error  has  the  decimal  on  the  wrong  side 
of  the  digit:  i.e.,  ,4->4.)).  The  peak  position 
reflects  the  total  stack  travel  time  while  the 
amplitude  decrease  reflects  the  loss.  The  5  Hz  peak 


may  be  due  to  the  spectral  contamination  by  the  surface 
waves  present  in  the  surface  records.  Further  analysis 
will  require  careful  windowing  to  isolate  the  shear 
wave.  The  positions  of  the  peaks  are  indicative  of  the 
travel  time  while  the  relative  peak  heights  reflect 
the  loss. 

Following  Figure  ?  all  nine  acceleration  records 
are  shown  along  with  their'  spectra. 

III.  CONCLUSION 

In  this  section  some  success  has  been  demonstrated 
in  modeling  the  response  of  a  soil  structure.  The 
method  used  an  extension  of  the  Haskell-Thompson 
matrix  method  presented  in  Section  1.  In  particular 
it  was  shown  that  a  significant  amount  of  attenuation 
was  necessary  in  order  to  compare  amplitudes.  It  was 
also  shown  that  a  large  portion  of  significant  motion 
was  due  to  lateral  propagation  and  therefore  not  present 
in  the  calculated  response. 

In  addition  the  observed  ratios  (surface  to  bedrock) 
showed  good  agreement  between  components.  Also,  it  was 
found  that  the  observed  ratios  compared  favorably  with 
the  calculated,  both  in  peak  position  and  in  peak 
duminition  with  frequency. 


ace.  hx  component 


I 


The  following  nine  figures  represent  the 
acceleration  records  arid  their  spectra  (dashed  line 
is  noise  spectra  windowed  prior  to  signal)  as  recorded 
by  the  Richmond  Fiel  Station  vertical  array  for  the 
Briones  Hills  event.  Components  are  labeled.  The 
poor  bedrock  (ZB)  record  was  the  result  of  telementry 
difficulty  which  has  subsequently  been  repaired. 
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CHAPTER  2 
SURFACE  WAVES 

I.  FORMULATION 

The  development  presented  in  Chapter  1  which 
incorporates  attenuation  in  an  exact  manner  into  a 
Haskell-Thompson  formulation  is  extended  here  to 
surface  waves.  The  extension  follows  quite  naturally 
for  the  period  equation  as  the  horizontal  wave  number 
now  becomes  complex  in  order  to  accommodate  the  surface 
wave  quality  factor.  This  now  requires  the  determin¬ 
ation  of  complex  roots  of  the  complex  determanental 
equation.  The  details  are  outlined  in  Appendices 
1  and  3* 

II.  FORWARD  PROBLEM 

In  this  section  the  surface  wave  formulation  (Love 
sind  Rayleigh  waves)  will  be  applied  to  a  soil  structure 
and  a  crust  and  upper  mantle  structure.  It  is  thought 
that  this  will  demonstrate  in  which  parts  of  the  earth 
and  over  which  period  ranges  anelastic  attenuation  may 
be  significant. 

A.  Soils 

A  typical  soil  structure  is  that  of  San  Francisco 
Bay  mud  and  is  shown  in  Table  1  of  Chapter  1.  Velocity 
and  density  values  were  measured  from  samples  taken 


from  a.  borehole  located  at  the  Richmond  Field  Station 
while  Q  values  represent  a  best  guess  (Chapter  1, 
Section  I). 

1.  Love  Waves 

Figure  1  shows  the  dispersion  and  attenuation 
parameters  for  Love  waves  over  the  soil  structure. 

The  phase  velocities  C  (E,  elastic;  L,  including  loss) 
are  very  similar  (no  causality  corrections  applied) 
and  shows  the  anelastic  phase  velocity  can  actually 
be  greater  than  the  elastic.  The  group  velocity  U£ 
is  the  elastic  because  the  Q  structure  (~10)  is  too 
low  to  apply  the  variational  method  of  Appendix  2 
meaningfully. 

The  Qt  curve  is  the  phase  quality  factor  and 
has  a  frequency  dependence  similar  to  the  elastic 
group  velocity  curve.  It  is  most  essential,  since  the 
group  quality  factor  is  the  physically  meaningful 
parameter,  to  develop  some  reliable  means  of  obtaining 
the  attenuating  group  velocity  for  highly  attenuating 
media . 

2.  Rayleigh  Waves 

In  Figure  2  are  shown  the  dispersion  and  attenu¬ 
ation  parameters  for  Rayleigh  waves  over  the  same  soil 
structure.  In  this  plot  we  note  some  interesting 
features.  The  two  phase  velocities,  elastic  and  loss, 
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are  asymptotic  at  the  longer  periods  but  diverge  sig¬ 
nificantly  toward  the  shorter  periods.  The  attenu¬ 
ating  phase  velocity  (C^)  is  significantly  greater 
than  the  elastic  at  the  shorter  periods  and  actually 
very  nearly  corresponds  to  the  first  overtone  elastic 
phase  velocity  for  periods  less  than  about  0.15 
seconds  (near  the  Airy  phase).  Meanwhile,  the  first 
overtone  attenuating  phase  velocity  corresponds  nearly 
to  the  elastic  fundamental  phase  velocity  over  this 
range.  This  suggests  that  a  mode  has  been  skipped 
in  the  calculation  procedure  but  this  does  not  appear 
to  be  the  case  because  the  associated  with  the  two 
modes  are  considerably  different  and  both  appear  to 
be  continuous  (a  much  finer  sample  interval  was  used 
in  this  range  to  check  continuity).  This  presents 
the  interesting  possibility  that  two  modes  may  possess 
the  same  phase  velocity  at  some  period  but  degeneracy 
is  avoided  through  distinct  attenuation  factors. 

This  may  be  the  result  of  the  rather  large  layer 
attenuation  significantly  perturbing  the  mode  shapes. 
Further  work  along  these  lines  is  needed.  These 
results  also  demonstrate  the  need  for  a  means  of 
accurately  calculating  the  group  velocity  (or  some 
more  meaningful  physical  parameter)  in  highly  attenu¬ 
ating  media. 

In  Figure  3  is  shown  the  phase  of  the  surface 
displacement  ratio  bo/VQ  for  the  elastic  and  loss 
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cases.  The  two  curves  are  significantly  different 
from  each  other,  particularly  in  the  region  near  the 
Airy  phase.  The  two  moduli  are  shown  in  Figure  4. 

It  is  interesting  that  the  magnitude  remains  nearly 
the  same  at  the  peak  while  the  loss  shifts  it  to 
shorter  periods  (towards  the  Airy  phase).  Clearly 
there  are  significant  differences  between  the  surface 
waves  for  elastic  and  highly  dissipative  media,  and  a 
complete  understanding  of  these  differences  will 
require  more  work. 

B.  Upper  Mantle 

The  upper  mantle  elastic  structure  which  is 
investigated  in  this  section  is  basically  that  shown  in 
Knopoff  and  Chang  (197?)  for  a  typical  oceanic  structure 
and  is  listed  in  Table  1.  An  oceanic  model  was  used 
because  it  is  appropriate  for  the  data  to  be  considered 
in  the  next  section,  which  was  collected  for  great 
circle  paths  that  were  over  70%  oceanic. 

The  attenuation  structures  were  chosen  to  be 
broadly  consistent  with  popular  models  yet  tailored 
to  give  a  chosen  starting  fit  for  the  inverse  problem 
considered  in  the  next  section.  They  are  therefore 
different  for  the  two  data  sets  considered  (see 
Section  III). 
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1.  Love  Waves 

Figure  5  shows  the  dispersion  and  attenuation 
parameters  for  Love  waves  over  the  upper  mantle  struc¬ 
ture.  Again  QL  is  the  phase  quality  factor  and  the 
group  quality  factor.  Both  the  phase  velocity  C  and 
group  velocity  U  include  attenuation  and  are  within 
0.01$  of  the  corresponding  elastic  results.  Group 
velocity  is  calculated  as  outlined  in  Appendix  2. 

No  dispersion  corrections  are  applied  (Kanamori  and 
Anderson,  1977)  since  the  effect  is  too  small  to  be 
noticable.  The  SH  earth  stretching  transformation 
was  applied  according  to  the  method  of  Schwab  and 
Knopoff  (1972). 

2.  Rayleigh  Waves 

Figure  6  shows  the  same  parameters  for  Rayleigh 
waves.  Again  dispersion  corrections  are  not  applied. 

The  P-SV  earth  stretching  transformation  was  applied 
according  to  the  method  of  Schwab  and  Knopoff  (1972). 

III.  INVERSE  PROBLEM 

In  order  to  apply  the  inversion  kernels  presented  in 
Appendix  I  on  observational  data,  an  initial  attempt 
was  made  at  inverting  400-25  sec.  fundamental  mode 
Love  and  Rayleigh  wave  data.  The  data  consist  of  phase 
velocity  and  phase  quality  factors  for  great  circle 
paths  through  Berkeley  recorded  on  a  broadband  system. 
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The  main  objective  of  this  exercise  is  not  to  derive 
a  realistic  upper  mantle  structure,  but  rather  to 
investigate  the  possible  difficulties  present  in  using 
the  exact  kernels  with  conventionally  extracted  amplitude 
data. 

A.  Data  Analysis 

For  the  Berkeley  data  the  following  procedure  was 
used.  The  record  tapes  are  filtered  to  remove  alias 
frequencies  and  then  digitized.  Seismograms  are  then 
group  velocity  windowed  between  4.50  and  4.30  km/sec 
for  Love  waves  and  between  4.00  and  3*45  km/sec  for 
Rayleigh  waves.  The  mean  is  removed  before  spectra 
are  taken.  Phase  velocity  is  then  determined  using 
the  smoothed  (Appendix  5)  phase  difference  for  two 
great  circle  paths  (Toksoz  and  Ben-Menahem,  1963). 

The  amplitude  spectra  is  smoothed  with  a  seven  point 
moving  average  and  then,  if  necessary,  further  smoothed 
by  eye.  The  phase  quality  factors  (defined  by  exp(-  ^q) 
C  is  the  phase  velocity)  are  determined  using  the  single 
station  method  (Kanamori,  1970). 

The  other  data  set  was  taken  from  Anderson  et  al . 
(1965)  with  both  the  Love  and  Rayleigh  wave  data 
augmented  at  the  short  and  long  periods.  For  Love  waves 
the  long  period  data  (T=400,360)  are  from  Press  et  al. 
(1961)  the  short  period  data  (T-39.55,  36.38,  24.92) 
are  from  two  sources;  the  velocity  data  are  synthetic 
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from  Case  122  of  Sykes  ejt  al.  (1962)  and  the  attenuation 
data  are  that  of  Tsai  and  Aki  (1969)*  The  Rayleigh  wave 
long  period  data  (T=400.00,  370.37)  are  from  Ben- 
Menahem  and  Toksoz  (1962).  The  short  period  data 
(T=39-29*  33.12,  24.63)  are  also  from  two  sources: 
the  velocities  are  synthetic  from  model  8099  of  Dorman 
et  al .  (I960)  while  the  attenuation  data  are  from 
Tsai  and  Aki  (1969) •  The  conversion  to  phase  quality 


factors  (Ql,  Qr)  has  been  made  according  to  the 

enuation  -  Q  =  0  U  group  velocity 

equation  Q  Qy  «L(R  c  phase  veiocity. 


B.  Inversion  Algorithm 

The  inversion  method  is  simply  a  damped  Gauss- 
Newton  algorithm  using  the  fundamental  decomposition 
theorem  (Lancsoz,  1961 ),  as  in  Appendix  1  with  synthetic 
data,  to  obtain  the  singular  values  and  the  inverse 
of  the  Jacobian.  An  overdetermined  system  is  used 
since  it  is  assumed  that  the  Q  starting  model,  although 
hopefully  close  enough  for  local  convergence,  probably 
is  in  considerable  error.  The  resolution  matrix 
(Jackson,  1972)  is  used  along  with  the  condition  number 
to  guide  in  structure  reduction  and  in  parameter 
constraint.  Standard  errors  of  the  parameter  adjust¬ 
ments  are  calculated  using  the  goodness  of  fit  as  an 
estimate  of  the  data  variance.  In  order  to  avoid  fixing 
parameters  at  the  wrong  value,  only  those  parameters 
whose  corrections  failed  to  stabilize  were  held  constant. 


This  allows  the  other  parameters  to  adjust  themselves 
somewhat  more  freely  without  compensating  for  the  fixed 
erroneous  value.  This,  of  course,  results  in  somewhat 
unreasonable  values  of  the  less  well  determined  parameters. 
However,  the  standard  errors  are  usually  large  thus 
keeping  the  parameters  within  acceptable  bounds. 

During  the  course  of  many  trial  and  error  inver¬ 
sions,  it  was  found  that  a  critical  factor  is  the  choice 
of  layering.  This  makes  a  very  good  case  for  either  of 
two  alternatives:  l)  the  use  of  the  above  procedure 
with  thickness  derivatives  to  aid  in  the  choice  of 
layering;  2)  the  underdetermined  system  coupled  with 
a  suitable  resolution  and  or  stabilization  method. 

In  favor  of  the  underdetermined  system  the  point  should 
be  made  that  perhaps  the  same  layering  may  not  be 
optimum  for  both  a  velocity  and  attenuation  inversion. 
Adjustments  for  this  could  easily  be  made  in  the  under¬ 
determined  system  simply  by  using  different  resolving 
kernels  (Knopoff  and  Jackson,  1972).  This  must,  of 
course,  be  weighed  against  the  heavy  dependence  upon 
starting  models  implicit  in  the  underdetermined  system. 

C.  Love  Waves 

Figure  7  shows  the  spectra  (a,  unsmoothed  and 
noise;  b,  smoothed)  of  the  G  pulses  3  through  6  for 
the  Tang  Shan  event  (7/27/76,  A =88°,  M=7.5).  For 


the  Ql  data  (Figure  9)  the  harmonic  mean  of  the  two 
ratios  6/4  and  5/3  was  used.  The  phase  data  (Figure  9) 
were  taken  from  the  6/4  ratio  only  as  it  seemed  the  best 
behaved.  The  phase  difference  (A/)  of  the  observed 
and  smoothed  is  shown  in  Figure  8. 

In  Table  1  is  shown  the  starting  model  used  in  the 
inversions  (Anderson  et  al.  augmented  data  has  a 
different  Q  starting  model  and  is  listed  in  appropriate 
tables).  The  model  (elastic  parameters)  is  a  typical 
oceanic  structure  (Knopoff  and  Chang,  1977)  slightly 
modified  to  aid  the  Love  wave  convergence  since  it  is 
a  poorer  resolvent  than  Rayleigh  waves.  The  Q 
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structure  was  designed  from  trial  and  error  to  give 
the  closest  starting  fit  while  being  as  generally 
consistent  as  possible  with  other  results.  The  notable 
exception  being  the  lack  of  a  low  Q  layer  coinciding 
with  the  low  velocity  layer.  Any  attempt  to  start  a 
low  Q  layer  there  resulted  in  either  divergence  or 
a  definite  increase  and  slow  convergence.  It  simply 
seems  to  be  incompatible  with  this  data.  The  crustal 
Q  was  fixed  at  500  while  the  half  space  Q  was  fixed  at 
200.  The  crustal  and  half  space  velocities  were  also 
fixed  due  to  the  limited  bandwidth.  Allowing  the  half 
space  parameters  to  vary  did  not  significantly  affect 
the  layer  above  but  resulted  in  less  reasonable  half 
space  parameters  (particularly  velocities). 


In  Table  2  are  shown  the  inversion  results  using 
the  Berkeley  data.  The  most  significant  result,  which 
appears  to  be  general,  is  the  incompatibility  between 


the  velocity  and  amplitude  data  with  respect  to  the 
velocity.  That  is,  inclusion  of  the  derivative  ^  ^L,R 

^trs 

generally  results  in  an  unstable  system.  The  ^ 
residuals  are  reduced  at  the  expense  of  the  C  residuals 
and  the  system  becomes  unstable  after  several  iterations 
(thirty-one  is  the  maximum  number  of  iterations  in  all 
cases).  A  look  at  sections  (B)  and  (D)  of  Table  2 
reveals  the  unfortunate  dilemma.  Section  (B)  shows  the 
results  without  the  derivative  ^ .  The  velocity 

J\T  s 

structure  is  reasonable  with  the  sum  square  residuals 
reduced  by  a  factor  of  ten  while  the  Q  structure  shows 
a  high  Q  lid  with  a  broad  low  Q  zone  below  the  low 
velocity  layer  and  the  sum  square  residuals  reduced  by 
1/3.  The  large  discrepancy  between  the  residual 
reduction  in  the  velocity  and  in  the  attenuation  might 
be  an  indication  that  the  amplitude  data  is  also  not 
representative  of  the  Q  structure.  Section  (D)  shows 
the  results  including  the  derivative  which 

(JYs 

became  unstable  after  the  ninth  iteration.  Note  the 
phase  velocity  residuals  have  not  been  reduced  while 
the  Q  residuals  are  reduced  by  another  1/3  over  section 
(B).  A  look  at  the  standard  errors  of  the  velocities 
(A  Vs)  for  both  cases  shows  generally  smaller  values 
when  the  derivative  ^ is  included  in  spite  of  the 


I 


much  larger  residuals.  This  indicates  more  information 
(orthogonal  to  )  as  discussed  in  Appendix  1; 

however,  a  look  at  the  two  velocity  structures  shows 
that  this  information  is  not  compatible.  In  section 
(C)  is  shown  the  same  conditions  as  in  (B)  without  the 
dispersion  correction  applied.  The  result  is  a  more 
pronounced  low  velocity  layer  (Hart  et  al .  ,  1976)  with 
generally  lower  velocities  throughout.  The  Q  structure 
is  also  somewhat  different,  as  expected,  but  still 
generally  indicates  a  broad  low  Q  zone  below  the  low 
velocity  layer.  Figure  10  shows  the  corrected  data 
along  with  the  phase  velocity  and  phase  quality  factor 
for  the  derived  model  of  section  (B). 

In  Table  3  are  shown  the  inversion  results  with 
the  augmented  Anderson  et  al.  Love  wave  data.  Here 
only  two  cases  are  shown:  (B)  without  the  derivative 

and;  (C)  including  the  derivative  with  dispersion 

t)Vs 

corrections  applied  in  both  cases.  The  starting  model 
(A)  has  a  different  Q  structure  which  was  based  on 

b 

the  Ql  data  (data  listed  in  Table  6).  The  general 
pattern  is  very  similar  to  the  Berkeley  data.  Neither 
the  velocity  nor  the  Q  structure  are  si/qiif icantly 

o 

different;  however,  thin  structure  is  generally 
lower  with  a  much  more  attenuating  lid.  Again,  there 
is  no  low  Q  zone  corresponding  to  the  low  velocity 
layer  but  rather  a  very  broad  low  Q  zone  below  the 
layer.  Also,  a  look  at  section  (C)  which  includes 
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the  derivative  ^  reveals  a  similar  pattern  to  the 

an 

Berkeley  data.  The  residuals  were  reduced  at  the 
expense  of  the  increased  C  residuals,  a  different 
velocity  structure,  and  an  unstable  system  after  six 
iterations.  This  indicates  that  the  augmented  Anderson 
et  al .  (1965)  Ql  data  is  also  inconsistent  with  respect 
to  velocity.  Also,  the  poor  QL  residuals  in  section  (B) 
indicate  that  the  Q-^  data  may  not  be  consistent  with 
respect  to  a  Q  structure.  The  data  and  fit  are  shown 
in  Figure  11. 

D.  Rayleigh  Waves 

Figure  12  shows  the  spectra  (a,  unsmoothed  and 
noise;  b,  smoothed)  for  R^  and  R^  for  the  Indonesian 
event  (8/19/77 ,  A=ll8° ,  M=8.0).  Due  to  noise  and  system 
problems  only  R^  and  R^  were  used  for  QR  while  R2  and 
were  used  for  the  phase  data.  Figure  13  shows  the 
phase  difference  A/  and  the  smoothed  fit. 

Table  4  shows  the  results  for  the  Berkeley  data. 

The  starting  model  (A)  is  the  same  as  in  the  Love  wave 
case  using  the  Berkeley  data.  The  starting  Q  structure 
is  simply  2.25  and  is  an  inversion  paramr+er  set 

o 

(Appendix  2).  The  V  structure  is  varied  by  maintaining 

X”' 

a  fixed  starting  Poisson's  ratio.  The  density  remains 
fixed.  For  the  Rayleigh  wave  inversions  it  was  found 
necessary  to  vary  the  half  space  parameters  to  allow 
the  layer  above  to  converge  to  a  reasonable  velocity. 
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Twenty-four  was  the  maximum  number  of  iterations 
used  in  all  Rayleigh  wave  inversions. 

In  section  (B)  of  Table  4  are  shown  the  results 
without  the  derivative  $  ^R.  The  residuals  have  been 

TTi 

reduced  nicely;  one  order  for  Qp  and  two  orders  for  C. 

However,  both  the  Q  and  V  structures  show  oscillations 

s  s 

which  are  probably  due  to  non-optimum  layering.  Since 
the  layering  was  guided  by  the  Love  wave  resolution 
matrix  only  (because  they  have  poorer  resolution  than 
Rayleigh  waves  (Appendix  1))  and  since  the  resolution 
matrix  is  made  up  of  partial  derivatives  which  have 
different  shapes  for  the  wave  types,  a  different  layering 
is  probably  required  for  the  case  of  Rayleigh  waves. 

This  also  makes  another  argument  for  the  use  of  an 
underdetermined  system.  Figure  14  shows  the  data  and 
the  model  fit. 

In  section  (C)  of  Table  4  are  shown  the  results 
including  the  derivative  d  ^R ■  Again,  the  results  are 

<J  vs 

very  similar  to  the  Love  wave  inversions.  The  system 
becomes  unstable,  the  C  residuals  are  not  significantly 
reduced,  and  the  solution  becomes  more  oscillatory  in 
the  velocity.  This  indicates  that  the  and  C  are 
competing  for  different  velocity  structures.  Again, 
since  the  QR  residual  in  section  (B)  has  not  been 
reduced  as  significantly  as  the  C  residual,  it  may 
indicate  that  the  data  is  not  representative  of  the 


attenuation  in  the  earth. 
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Turning  now  to  the  augmented  Anderson  et  al.  (1965) 
data,  the  results  are  shown  in  Table  5«  Section  (A) 
shows  the  starting  model  which  is  identical  to  the  Love 
wave  case  for  the  same  data.  As  in  the  previous  case 
the  solution  is  oscillatory  though  stable  when  the 
derivative  t)  is  not  included  (section  (B))  and 
unstable  with  respect  to  velocity  when  the  derivative 
is  included  (section  (C)).  No  further  conclusions  should 
be  drawn  from  the  Rayleigh  wave  analysis,  especially 
with  regard  to  parameter  values. 

IV.  CONCLUSION 

Concerning  the  forward  calculations,  some  inter¬ 
esting  features  were  noted  for  Rayleigh  waves  in  high 
loss  material.  It  appears  that,  near  the  Airy  phase, 
the  loss  can  significantly  perturb  mode  shapes  such 
that  the  phase  velocity  can  become  multi-valued. 

Degeneracy  is  avoided  through  distinct  attenuation 
factors  for  each  mode.  These  observations  are  pre¬ 
liminary  and  clearly  more  work  needs  to  be  done.  Also, 
a  reliable  means  of  determining  group  velocity  for  such 
high  loss  material  must  be  developed. 

In  the  inverse  problem  investigations,  probably 
the  most  salient  result  is  the  incompatibility  of  the 
Q  data  with  the  C  data  with  respect  to  the  velocity 
structure.  The  incompatibility  is  expressed  through 


the  derivative  u. *L,R,  its  inclusion  rendering  the 
system  unstable  in  velocity  and  adversely  affecting 
the  C  residuals.  The  fact  that  the  QL  R  residuals 
are  much  larger  than  the  C  residuals  when  the  derivati  e 

j ^L,R  is  neglected  also  indicates  that  conventionally 

£?  \Ts 

extracted  surface  wave  Q  data  probably  are  in  consider¬ 
able  error.  The  error  is  probably  due  to  multi-path 
propagation  and  affects  the  spectral  modulus  to  a 
greater  degree  than  the  phase  (Pilant  and  Knopoff,  1964). 
A  method  of  accurate  amplitude  determination  which  seems 
hopeful  is  that  of  phase-matched  filters  (Herrin  and 
Goforth,  1977)  and  future  work  will  be  undertaken  along 
these  lines. 

Another  result  which  emerges  from  the  inversion 
investigations  is  the  possible  inadequacy  of  the  over¬ 
determined  system.  It  appears  that  it  would  be  most 
difficult  to  derive  a  layering  that  is  simultaneously 
optimum  for  velocity  and  attenuation  with  both  Love 
and  Rayleigh  waves.  To  optimize  the  Love  and  Rayleigh 
layering,  a  simultaneous  inversion  seems  appropriate, 
but  this  does  not  resolve  the  velocity  and  attenuation 
layering  problem.  Some  investigations  are  needed  into 
the  development  of  a  suitable  underdetermined  formu¬ 
lation  with  its  heavy  dependence  on  starting  models. 

Another  approach,  if  a  model  rather  than  local 
averages  is  desired,  might  be  the  use  of  a  suitable 
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non-linear  algorithm  such  as  the  Marquardt-Levenberg 
(Marquardt,  1963).  Perhaps  this  algorithm  would 
tolerate  a  larger  condition  number  (more  layers)  thereby 
reducing  the  layering  incompatibility. 
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CHAPTER  3 
SUMMARY 

Lately,  albeit  belatedly  (Jeffreys,  1965),  it 
has  become  apparent  that  the  earth  can  no  longer  be 
treated  as  an  elastic  body.  An  essential  point  of  this 
work  has  been  to  establish  that  it  is  equally  insuf¬ 
ficient  to  treat  it  as  a  perturbation  on  an  elastic 
body. 

An  exact  plane  layer  propagator  matrix  which 
includes  attenuation  was  presented.  Application  of 
the  formulation  to  continue  an  incident  elastic  wave 
through  a  realistic  soil  structure  to  the  earth’s 
surface,  demonstrated  the  importance  of  energy 
absorption  in  predicting  ground  motion.  The  minor 
importance  of  converted  waves  was  demonstrated  and 
it  was  further  shown  that  lateral,  propagation,  not 
included  in  the  formulation,  can  affect  the  duration 
of  motion.  In  addition  to  time  domain  analysis, 
spectral  ratios  (surface  to  bedrock)  show  fair  agree¬ 
ment  with  model  calculations. 

In  applying  the  matrix  formulation  to  surface 
wave  eigenvalue  and  surface  displacement  calculations, 
several  results  were  discussed.  It  appears  that 
attenuation  perturbs  Rayleigh  waves  significantly  more 
than  Love  waves  in  a  high  loss  structure.  In  par¬ 
ticular,  different  mode  phase  velocities  may  cross 


/ 


67 

near  the  Airy  phase.  Degeneracy  is  avoided,  however, 
through  distinct  quality  factors. 

Surface  wave  inversions  for  both  velocity  and 
attenuation  structures  using  an  overdetermined  system 
revealed  an  incompatibility  with  respect  to  the 
velocity  structure  between  phase  and  attenuation  data. 

)  Q 

In  particular  inclusion  of  the  derivative  v  L.R 

2  \r  s 

with  real  data  greatly  degraded  the  velocity  solution 
compared  to  inversions  using  only  the  derivative 

y  C  . .  Poor  attenuation  solutions  compared  with 

t)\fs 

velocity  solutions  further  indicates  the  inadequacy 
of  presently  used  attenuation  data.  Clearly  a  reliable 
formulation  is  needed  in  order  to  extract  meaningful 
amplitude  data  from  surface  wave  observations. 

Of  particular  concern  in  dealing  with  an  over¬ 
determined  system  which  also  bears  on  the  quality  of 
attenuation  solutions,  is  the  choice  of  layering. 

This  appears  to  be  critical  at  least  for  the  velocity 
structure.  The  same  layering  may  not  be  optimum  for 
both  Love  and  Rayleigh  waves  in  addition  to  both 
velocity  and  attenuation  structures. 

The  general  conclusion  of  this  thesis  is  that 
amplitude  information  is  most  essential  in  furthering 
our  understanding  of  both  the  velocity  and  attenuation 
structures  of  the  earth.  It  has  also  been  demonstrated 
that  realistic  calculations  with  earth  materials  require 
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TABLE  1 


Starting  model  for  both  Love  and  Rayleigh  wave 
inversions.  The  model  is  based  on  a  typical  oceanic 
structure  (Knopoff  and  Chang,  197?).  The  low  velocity 
layer  thickness  was  increased  20  km  at  the  expense  of 


the  layer 

below  to  aid  convergence 

in  the 

Love 

wave 

inversions 

(Berkeley  data) 

thickness 

Vp(k»/sec) 

Vg(km/sec) 

(cgs) 

Qp 

Qs 

(km) 

2.10 

1.00 

2.10 

1200 

500 

1.0 

6.41 

3-70 

3-07 

1200 

500 

5.0 

8.10 

4.65 

3.40 

4500 

2000 

70.0 

7.60 

4.15 

3.40 

425 

200 

130.0 

8.80 

^.75 

3.65 

275 

125 

240.0 

9.80 

5.30 

3.98 

225 

100 

200.0 

11.15 

6.20 

4.43 

225 

100 

400.0 

11.78 

6.48 

4.63 

350 

150 

240.0 

12.02 

7.20 

4.71 

450 

200 

cs«» 

Fluid  layer  included  in  Rayleigh  calculation 
1.50  0.00  1.00  1200 


4.0 
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TABLE  2 


Results  of  Love  wave  inversions  with  Berkeley  data: 
A)  starting  models  B)  inversion  without  the  derivative 


with  dispersion  correction  applied;  C)  inversion 

«L 


i  Vs 

without  both  the  derivative  -1-  *  and  dispersion 

'J  V  S 

correction;  D)  inversion  including  both  the  derivative 
Qt 

-r-=-  and  dispersion  correction.  (AQ_,  A  V_  =  standard 
» s  s  s 

error;  SSRQ,  C  »  sum  square  residual  Qr  l’  CR  L!  ^ ’ 
(F)  =  initial,  final;  NI  =  number  of  iterations). 


Layer 

^s 

AQ 

s 

*Vs 

SSRQ 

SSRC 

1 

500 

1.00 

2 

500 

3-70 

3 

2000 

4.65 

4 

200 

^.15 

A 

5 

125 

4.75 

6 

100 

5-30 

7 

100 

6.20 

8 

150 

6.48 

9 

200 

7.20 

1 

500 

1.00 

2 

500 

3.70 

3 

2289 

20050 

4.47 

0.02 

4 

165 

129 

4.32 

0.03 

(I) 

6.0 

7.0x10 

B 

5 

144 

100 

4.62 

0.02 

(F) 

2.0 

4.0x10 

6 

84 

142 

5.37 

0.08 

7 

80 

159 

6.77 

0.10 

8 

126 

85 

7.02 

0.07 

9 

200 

7.20 

1 

500 

1.00 

2 

500 

3.70 

3 

2090 

47413 

4.47 

0.06 

4 

196 

363 

4.29 

0.06 

C 

5 

124 

46 

^■55 

0.02 

(I) 

6.0 

2.0x10 

6 

102 

156 

5-36 

0.07 

(F) 

2.0 

5.0x10 

7 

87 

167 

6.57 

0.10 

8 

144 

82 

b .  80 

0.06 

9 

200 

7.20 

/ 


TABLE  2  (continued) 


Layer 

aQ3 

V 

_ a 

SSRQ 

SSRC 

1 

500 

1.00 

2 

500 

3.70 

3 

194 

5854 

4.42 

0.00 

4 

198 

12? 

4.44 

0.00 

(I) 

6.0 

7.0x10 

D  5 

122 

51 

4.46 

0.01 

(F) 

1.0 

7.0x10 

6 

97 

56 

5-70 

0.08 

7 

99 

22 

6.23 

0.10 

3 

151 

13 

6.52 

0.08 

9 

200 

7.20 
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TABLE  3 


Results  of  Love  wave  inversions  with  augmented 
Anderson  et  al.  (1963)  data:  A)  starting  model 
(different  Qg  structure  only);  B)  inversion  without 


t)  ®L 


the  derivative  *»  ***;  C)  inversion  including  the 

„  .  isJ7* 

derivative 

o  fs 

all  cases. 


Dispersion  corrections  applied  in 


Layer 

^s 

*Qs 

la 

SSRQ 

SSRC 

1 

500 

1.00 

2 

500 

3.70 

3 

500 

4.65 

4 

100 

4.15 

A  5 

100 

4.75 

6 

100 

5-30 

7 

100 

6.20 

8 

100 

6.48 

9 

200 

7.20 

1 

500 

1.00 

2 

500 

3.70 

3 

423 

3470 

4.99 

0.06 

4 

126 

77 

4.18 

0.01 

B  5 

88 

124 

4.43 

0.04 

(I) 

14 

3-0x10 

6 

108 

393 

5 . 66 

0.11 

(F) 

12 

4.0x10 

7 

88 

429 

6.18 

0.19 

3 

129 

203 

6.33 

0.10 

9 

200 

7.20 

1 

500 

1.00 

2 

500 

3.70 

3 

505 

871 

4.66 

0.02 

4 

102 

113 

4.37 

0.00 

C  5 

101 

132 

4.56 

0.03 

(I) 

14 

3.0x10 

6 

104 

76 

5.51 

0.14 

(F) 

7 

2.0x10 

7 

104 

40 

5-94 

0.13 

8 

156 

21 

6.21 

0.14 

9 

200 

94 

7.20 

TABLE  4 


Results  of  Rayleigh  wave  inversions  with  Berkeley 
data:  A)  starting  model,  B)  inversion  without  the 
derivative  -^~5;  C)  inversion  including  the  derivative 
Jv"""s'  DisPersion  corrections  applied  in  all  oases. 


1  500 

2  500 

3  500 

4  2000 

A  5  200 

6  125 

7  100 

3  ioo 

9  150 

10  200 


V 

-JB 


1200 

1200 

1200 

4500 

425 

275 

225 

225 

350 

450 


0.00 
1.00 
3.  ?0 

4.65 

4.15 

4.75 

5.30 

6.20 

6.43 

7.20 


1.52 

2.10 

6.41 

8.10 

7.60 

8.30 

9.80 

11.15 

II.78 

12.82 


SSSQ  SSRC  NI 


1 

2 

3 

4 

3  5 
6 

7 

8 
9 

10 


500 

5C0 

<00 

366 

172 

69 

395 


224 


27 

67 


1  500 

2  500 

3  500 

4  745 

C  5  126 

6  77 

7  230 

0  176 

9  70 

10  194 


1328 

18956 

0.00 

1329  152452 

1.00 

1304 

40199 

3.70 

104 

4604 

314 

4.21 

6 

394 

79 

4.08 

5 

322 

35 

4.99 

436 

197 

24 

4.66 

37 

66 

3 

6.73 

10 

229 

1 

7.32 

1 

447 

0 

7.35 

1.52 
2.10 
6.41 
0.01  7.33 

0.01  7.49  (I)  4.0 

0.02  9.24  If)  .4 

0.05  8.62 

0.06  12.00 
0.08  13.31 

0.00  13.09 


1093  5225  0.00 

1095  38155  1.00 
1161 
217  4409 
8  538 


13  322 

194  229 

23  217 

18  341 

1  449 


7962  3.70 

324  4.39 

62  3.80 


.  M-  ( 

4  6.36 

14  7.60 

0  3.28 


1.52 
2.10 
6.41 
0.06  7.64 

0.03  6.96 

O.C7  10.62 
0.13  8.26 

0.24  11.44 

0.33  13.82 

0.10  14.74 


(I)  4.0 

(F)  .6 


7.0x10”? 
5.0x10”*  24 


7.0xJ0'2  24 

2.0x10”^ 


i 
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TABLE  5 

Results  of  Rayleigh  wave  inversions  with  augmented 

Anderson  et  al.  (1965)  data:  A)  starting  model  (sane 

as  Love  wave  in  Table  4);  B)  inversion  without  the 
v  ®R 

derivative  j  C)  inversion  including  the  derivative 
p*  «r  .  3 

-y-~.  Dispersion  corrections  applied  in  all  cases. 


Layer 

1 

2 

3 

4 

A  5 
6 
7 
3 
9 

10 

1 

2 

3 

4 

3  5 
6 
7 
3 
9 

10 

1 

2 

3 

4 

C  5 
6 

7 

8 
o 

10 


% 

^0 

!s 

!  !a 

500 

1200 

o.co 

1.52 

500 

1200 

1.00 

2.10 

500 

1200 

3.70 

6.41 

500 

1200 

4.65 

8.10 

100 

225 

4.15 

7.60 

100 

225 

4.75 

3.80 

100 

225 

5.30 

9-80 

100 

225 

6.20 

11.15 

150 

350 

6.48 

11.78 

200 

450 

7.20 

12.82 

500 

1348 

39573 

0.00 

1-52 

500 

1341 

144248 

1.00 

2.10 

500 

982 

67967 

3.70 

6.41 

2017 

44 2  6 

1040 

790  9 

4.27 

0.03 

7.71 

91 

24 

363 

1082 

4.66 

0.03 

8.53 

213 

131 

309 

612 

4.28 

0.04 

7.94 

118 

112 

274 

190 

6.35 

O.39 

11.73 

503 

2106 

141 

221 

5.48 

0.31 

9- 36 

103 

79 

242 

ne 

7.63 

1.47 

13.88 

87 

60 

431 

3 

8.33 

0.14 

14.33 

500 

1175 

94 

0.00 

1.52 

500 

1133 

663 

1.00 

2.10 

500 

1122 

254 

3.70 

6.41 

165 

14 

1339 

14 

4.23 

0.09 

7-36 

34 

6 

296 

9 

4,74 

0.01 

8.69 

196 

15 

366 

3 

4.14 

0.06 

7.66 

149 

11 

364 

3 

6.28 

O.23 

11.62 

325 

3 

63 

1 

5-75 

0.16 

10.33 

112 

6 

211 

1 

5-73 

0.36 

10.51 

63 

0 

441 

0 

3.15 

0.02 

14.52 

(I) 

(F) 


SSRQ  3 SRC  NI 


8.0  2 .0x10?  24 

2.0  Z.OxlO'-5 


8.0  2 . 0x10 ",  24 

.1  3 


/ 
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TABLE  6 


Data  used  in  inversions:  A)  Berkeley  data: 
B)  augmented  Anderson  et  al.  (1966)  data. 


A) 


Love 


Rayleigh 


T 

«L 

C 

C^x 

409.6 

88 

5-609 

5.695 

341.3 

99 

5.349 

5-417 

292.6 

111 

5-177 

5.234 

256.0 

120 

5.055 

5.105 

227.6 

128 

4.964 

5.008 

204.8 

135 

4.894 

4.934 

186.2 

141 

4.838 

4.875 

170.7 

145 

4.792 

4.827 

157.5 

150 

4.755 

4.787 

146.3 

156 

4.723 

**•753 

136.5 

161 

4.696 

4.725 

128.0 

166 

4.674 

4.701 

120.5 

169 

4.654 

4.679 

113.8 

174 

4.637 

4.662 

107.8 

177 

4.622 

4.646 

102.4 

181 

4.609 

4.632 

97.5 

185 

4.598 

4.619 

93-1 

189 

4.588  ' 

4.609 

89.0 

192 

4.579 

4.599 

85-3 

196 

4.572 

4.592 

81.9 

198 

4.565 

4.584 

78.8 

201 

^•559 

4.577 

75-3 

203 

<*•553 

4.571 

73.1 

206 

4 , 548 

4.565 

70.6 

208 

4.543 

4.560 

68.3 

209 

4.538 

^.555 

66.1 

211 

4.534 

4.544 

64.0 

213 

4.529 

4.545 

62.1 

216 

^.525 

4.540 

60.2 

217 

4.520 

4.535 

T 

qr 

C 

384.0 

96 

5.895 

5.976 

341.3 

108 

5-504 

5-569 

307.2 

116 

5-229 

5.285 

279-3 

124 

5.026 

5.075 

256.0 

129 

4.869 

4.914 

236.3 

131 

4.746 

4.788 

219.4 

131 

4.646 

4.686 

204.3 

130 

4.563 

4.602 

192.0 

128 

4.493 

4.531 

j.30 . 7 

124 

4.433 

4.472 

170.7 

120 

4.382 

4.420 

161.7 

117 

4.336 

4.374 

153.6 

114 

4.296 

4.334 

146.3 

112 

4.259 

4.297 

139.6 

111 

4.225 

4.261 

133.6 

110 

4.194 

4.231 

128.0 

110 

4. 165 

4.202 

122.9 

111 

4.138 

4.174 

118.2 

113 

4.112 

4.146 

113.8 

116 

4.088 

4.120 

109.7 

120 

4.064 

4.095 

105.9 

124 

4.042 

4.071 

102.4 

130 

4.020 

4.048 

99.1 

136 

3.999 

4.025 

96.0 

143 

3.979 

4.004 

93.1 

149 

3.961 

3.983 

90.4 

156 

3.942 

3-964 

87.8 

164 

3.925 

3.945 

85.3 

172 

3.909 

3.928 

83.O 

180 

3.894 

3.912 
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TABLE  6  (continued) 


58.5  222 

56.9  22? 

55.4  233 

53.9  23? 

52.5  244 

51.2  250 

49.9  259 

48.8  2?1 

47.6  281 

46.6  294 

45.5  311 

44.3  329 

43.6  350 

42.?  36? 

41.8  375 

40.9  380 

40.2  386 

34.7  466 

30.1  49? 

26.6  547 


4.516 

4.530 

80.8 

18? 

4.511 

4.525 

78.? 

193 

4.507 

4.521 

76.8 

200 

4.503 

4.516 

71.6 

220 

4.499 

4.511 

68.3 

234 

4.495 

4.50? 

65 .4 

243 

4.491 

4.503 

62.? 

253 

4.488 

4,499 

59.1 

266 

4.485 

4.495 

56.9 

270 

4.482 

4.492 

53.9 

280 

4.497 

4.489 

51.2 

282 

4.478 

4.486 

48.8 

275 

4.476 

4.484 

66,5 

280 

4.475 

4.482 

44 .5 

283 

4 ,4?4 

4.481 

4  2.? 

288 

4.473 

4.480 

40.4 

279 

4.4?2 

4.479 

36.1 

320 

4.465 

4.470 

29.3 

321 

4.452 

4.457 

25.2 

369 

4.455 

4.459 

C 


3.880 

3.86? 

3.853 

3.82? 

3.814 

3.306 

3.802 

3.801 

3.803 

3.808 

3.812 

3.814 

3.812 

3.806 

3.799 

3.789 

3.782 

3.823 

3.73? 


3.987 

3.883 

3.8?1 

3.841 

3.827 

3.313 

3.813 

3.811 

3.813 

3.817 

3.821 

3.823 

3.820 

3.815 

3.80? 

3.79? 

3.788 

3.788 

3.73? 


B) 


400.0 

360.0 

333-3 

312.5 

294.1 
277.8 

263.2 


250.0 

238.1 

227.3 


112  5.500 

110  5.380 

113  5.307 

114  5.243 

115  5.185 

113  5.134 

112  5.088 

U3  5-046 

114  5.008 

116  4.972 


5.565 

5-443 

5.366 

5-300 

5-240 

5.189 

5.142 

5.098 

5.05? 

5.021 


^00.0  138 

370.4  147 

333.3  189 

312.5  170 

294.1  156 

277.8  157 

263.2  149 

250.2  142 

238.1  144 

227.3  138 


5.955  6.042 
5.760  5.811 
5.548  5. 58I 
5.383  5.422 
5.251  5.292 
5.132  5.171 
5.012  5.052 
4.91?  4.958 
4.824  4.863 
4.754  4,793 


217.4  115 

208.3  115 
200.0  116 

192.3  117 

185-2  11? 

178.6  118 

172.4  118 

166.7  119 

I6I.3  116 

156.2  117 


4 . 940  4 . 988 
4.911  d.959 
4.885  4.931 
4.861  4.906 
4.840  4.885 
4.821  4.865 
4.805  4.848 
4.790  4,332 
4.776  4.519 
4.761  4.303 


217.4  140 

208.3  136 

200.0  13? 

192.3  136 

185.2  136 

178.6  137 

172.4  133 

166.7  I38 

161.3  140 

156.2  141 


4.674  4.712 
4.616  4.656 
4.569  4.606 
4.51?  4.554 
4.469  4.504 
4.436  4.470 

4,40?  4.441 
4.374  4.4o? 
4.348  4.3S0 
4.319  4. 3  so 


/ 
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TABLE  6  (continued) 


T 

C 

151.5 

116 

4.748 

4.789 

125.0 

101 

4.690 

4.734 

113.6 

102 

4.673 

4.715 

104.2 

106 

4.662 

4.701 

96.2 

113 

4.652 

4.688 

39.3 

120 

4.650 

4.683 

83.3 

122 

4.627 

4.658 

78.1 

118 

4.613 

4.644 

73.5 

106 

4.603 

4.637 

69.4 

100 

4.595 

4.630 

65.8 

93 

4.538 

4.624 

62.5 

95 

4.577 

4.612 

59.5 

95 

4.558 

4.592 

56.8 

96 

4.539 

4.572 

54.3 

98 

4.522 

4.554 

52.1 

101 

4.519 

4.549 

50.0 

105 

4.49s 

4.526 

39.6 

164 

4.490 

4.506 

36.4 

246 

4.480 

4.490 

24.9 

492 

4.460 

4 .464 

T 

qr 

C 

CUrx 

151.5 

142 

4.299 

4.329 

125.0 

132 

4.250 

4.281 

113.6 

140 

4.196 

4.224 

104.2 

145 

4.173 

4.199 

96.2 

139 

4.167 

4.193 

89.3 

135 

4.159 

4.185 

83.3 

129 

4. 167 

4.194 

78.1 

128 

4.169 

4.195 

73.5 

128 

4.160 

4.185 

69.4 

127 

4.170 

4.195 

65.8 

125 

4. 163 

4.188 

62.5 

121 

4.159 

4.184 

59.5 

117 

4.155 

4.180 

56.8 

115 

4.150 

4.175 

5^-3 

116 

4.145 

4.169 

51.2 

118 

4.141 

4.164 

50.0 

122 

4.136 

4.158 

39.3 

151 

3.985 

4.000 

33-1 

170 

3.988 

4.000 

24.6 

166 

3.980 

3.990 
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40.9 


22.8 


120  240  360  SEC. 

Plots  of  the  smoothed  G  pulse  spec t rums .  Data  shown  to  40  sec 
only  as  three  shorter  period  data  were  later  added. 
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Fig.  12  a. 


Plots  of  the  Rayleigh  wave  spectrums.  Dashed 
curve  is  a  noise  sample  windowed  leading  the 
wavetrain  sampled. 
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APPENDIX  1 

INVERSION  OF  SURFACE  WAVE  DATA  FOR  VELOCITY  AND 
ANELASTICITY  USING  EXACT  KERNELS 


106 

I.  INVERSION  OF  LOVE  WAVE  DATA  FOR  VELOCITY 
AND  ANELASTICITY  USING  EXACT  KERNELS 

by 

Walter  Silva 

Department  of  Geology  and  Geophysics 
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ABSTRACT 

The  general  problem  of  inverting  Love  wave 

dispersion  and  amplitude  data  to  obtain  a  velocity 

and  Q  structure  is  considered.  A  formulation  is  used 
s 

which  incorporates  attenuation  into  the  Haskell -Thompson 

matrix  method  in  an  exact  manner  and  thus  retains  the 

inherent  non-linearity  in  the  anelasticity .  The 

resulting  exact  inversion  kernels  allow  simultaneous 

inversion  for  velocity  and  intrinsic  attenuation 

parameters.  The  method  is  applied  to  synthetic  data 

which  allows  a  comparison  to  be  made  with  inexact 

kernels.  The  results  indicate  that  the  use  of  inexact 

kernels  may  introduce  spurious  oscillations  into  the 

Q  structure  and  that  a  simultaneous  inversion  can  be 
s 

more  stable  than  inverting  for  velocity  alone. 


10? 
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INTRODUCTION 

Surface  waves  provide  an  invaluable  tool  in  the 
study  of  the  earth's  interior.  Their  use  has  mainly 
been  in  dispersion  studies  to  infer  such  elastic  para¬ 
meters  as  velocity  and  density  in  structures  ranging 
from  soils  to  the  lower  mantle.  However,  in  recent 
years  there  has  been  an  increasing  interest  in  the 
anelastic  parameters  of  the  earth.  This  type  of  in¬ 
vestigation  can  supply  valuable  information  concerning 
material  properties,  structures,  and  temperature  dis¬ 
tributions  in  the  earth.  Of  most  recent  interest  is 
the  frequency  dependence  and  coupling  between  intrinsic 
velocity  and  attenuation.  With  this  in  mind,  it  now 
becomes  important  to  have  more  exact  methods  of  in¬ 
verting  the  data  for  both  velocity  and  attenuation  if 
any  meaningful  interpretation  is  to  result  from  such 
investigations . 

A  formulation  is  presented  here  using  an  exact 
generalization  of  the  Haskell-Thompson  matrix  method 
(Haskell,  1953 i  Silva,  1976b)  to  include  anelasticity 
in  inverting  Love  wave  dispersion  and  attenuation  data 
for  both  layer  velocity  and  attenuation.  Synthetic 
data  are  used  to  allow  a  comparison  to  be  made  between 
the  exact  formulation  and  the  approximate  linear 
theory  of  Anderson  and  Archambeau  (1964).  A  somewhat 
unstable  problem  is  considered  in  order  to  demonstrate 
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the  differences  and  to  represent  more  accurately  actual 
inversions,  since  surface  wave  attenuation  data  are 
typically  sparse  and  of  limited  accuracy  and  bandwidth. 

FORMULATION 

In  a  recent  paper  (Silva,  1976b)  a  theory  was 
presented  for  introducing,  in  an  exact  manner,  anelastic 
attenuation  into  a  Haskell-Thompson  formulation.  This 
is  applied  here  to  Love  waves.  The  earlier  treatment 
considered  plane  P  and  SV  waves  propagated  in  a  layered 
linear  viscoelastic  half-space.  In  particular, 
attenuating  (or  elastic)  layers  over  an  elastic  half¬ 
space  were  considered.  However,  due  to  the  boundary 
conditions,  the  component  of  attenuation  parallel  to 
the  interface  Ax  (which  describes  the  spatial  decay  of 
the  surface  wave)  is  constrained  to  be  continuous  along 
with  the  horizontal  wave  number  Px>  This  means  that 
each  layer,  along  with  the  half-space,  must  be  attenu¬ 
ating. 

Define  P  and  A  as 

X  X 


with  c  the  Love  wave  phase  velocity  and  where  Ax  has 
been  defined  in  this  form  to  accommodate  a  Love  wave 
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phase  quality  factor  (Brune,  1962).  The  attenu¬ 
ation  then  has  introduced  another  eigenvalue  'to 

be  determined  as  a  function  of  period.  In  the  appli¬ 
cation  to  surface  waves,  equation  (18)  of  Silva  (1976b) 
(considering  now  SH  only)  for  the  complex  vertical 
component  K of  the  complex  wave  number  K  must  be 
modified  to  the  following: 

Kz  =  [kZ  -  KxZp  R  £k?>R  £kx?  (2) 

=  _i  [KxZ  -  kZ]  *  R  fK?  <  R  fKx? 


where 


with  V  and  Q  the  homogeneous  shear  wave  velocity  and 
specific  attenuation  factor,  respectively.  The  charac¬ 
teristic  equation  is  then  obtained  in  the  usual  way 
(Haskell,  1953)  and  Kx  found  for  each  period  using 
Mueller's  method  (Conte  and  Debor,  1972).  In  Figure  1 
are  shown  C,  Q^,  and  U,  the  group  velocity,  for  the 
crustal  model  of  Table  1  which  represents  the  error 
free  data.  Group  velocity  is  calculated  using  a 


variational  method  applied  to  dissipative  media 
(Silva,  1976a). 
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Since  the  development  above  represents  an  exact 
forward  solution,  the  non-linearity  of  both  the  shear 
velocity  and  intrinsic  attenuation  is  retained.  In 
order  then  to  invert  the  phase  velocity  and  attenuation 
data  for  V  and  Q  a  Taylor  series  approximation  is 

0  S 

made,  keeping  only  the  linear  terms,  and  iterations 
performed.  We  then  have: 


9  c'X 


Sq, 


PVS 

QL 


A  vo  +  4 


i>v3  A  Vs  +  ?'« 


A  Q. 


A  Q< 


(3) 


where  £  c  and  & represent  the  difference  between  the 
observed  phase  velocity  and  Love  wave  attenuation 
parameters  at  each  period  and  those  calculated  from 
an  appropriately  close  starting  model.  The  derivatives 
are  calculated  in  the  following  way:  the  first 
is  an  analytical  derivative  using  an  extension  of 
Rayleigh's  principle  to  dissipative  media  (Silva, 
1976a),  while  the  remaining  are  calculated  numerically, 
using  a  10%  parameter  change.  Plots  of  the  derivative 
versus  period  for  each  layer  are  shown  in  Figures  2 
through  5*  The  plots  contain  a  number  of  interesting 
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features.  Figure  2  shows  the  derivative 


L 


be  compared  with  p 
relation, 


V 


c_ 

V. 


Qt 
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and  should 


(Figure  5)  which  confirms  the 


()  QI 

P  Q£ 


(4) 


which  is  derived  as  equation  (A2)  in  Appendix  4.  This 
is  also  the  relation,  applied  rather  to  inverse  Q, 
used  in  the  linearized  theory  (Anderson  and  Archambeau, 

1964) .  The  most  revealing  plots  are  those  of  «  -r— 

!)  Qt  ;QS 

(Figure  3)  and  -c ~rr~  (Figure  4)  which  confirm  equation 
»  vs 

(A3)  of  the  Appendix  4: 


h. 


Q  ^ 

s 


4  Qs2  Ql2 


P  Vs 


(5) 


In  Appendix  4  it  is  argued  that  in  equations  such  as 

Ql 


(3)  the  term  involving 
.  ;c 

involving  -v  rjr  and 


is  comparable  to  those 
This  demonstrates  the 


7TS  —  ;  qs  • 

large  dependence  of  surface  wave  amplitudes  on  the 

velocity  structure  in  plane  layer  models  and  strongly 

«?qL 

suggests  that  the  derivative  v  will  play  an 

V  *  g 

important  part  in  any  inversion  scheme.  Also  note  the 
lack  of  similarity  of  the  curves  in  Figures  3  and  4 
with  those  in  Figures  2  and  5  which  indicates  that 
they  contain  a  different  distribution  of  information 
and  thus  their  inclusion  should  add  to  the  stability 
of  the  inversion  process. 
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TESTS  WITH  SYNTHETIC  DATA 


The  crustal  model  shown  in  Table  1  was  used  as 
a  layer  over  half-space  in  a  previous  work  (Harkrider, 
1968)  to  demonstrate  the  results  of  a  variational 
formulation  for  surface  waves.  The  40  km  thick  crust 
was  subdivided  here  into  four  identical  10  km  thick 
layers  in  order  to  allow  the  inversion  process  greater 
freedom  in  parameter  adjustment.  Although  this  may 
seem  excessive  in  order  to  exaggerate  the  effect,  one 
should  realize  that  the  data  here  are  error  free  and 
the  real  model  is  known.  The  period  range  is  120  to  5 
seconds  with  5  second  intervals  and  the  non-geometrical 
dispersion  due  to  dissipation  (Burton,  1977)  is 
neglected . 

In  order  to  compare  the  linearized  inversion 
theory  (Anderson  and  Archambeau,  1964)  with  the  exact 
formulation,  both  inversions  were  carried  out  on  the 
same  data  (Figure  1)  generated  by  the  exact  theory 
using  the  crustal  model  of  Table  1.  Both  inversions 
used  the  method  of  singular  value  decomposition 
(Lancsoz,  1961)  with  the  non-linear  appropriately 
scaled  (Wiggins,  1972)  in  order  to  produce  a  consistent 
comparison.  The  results  of  the  inversions  are  shown 
in  Table  2. 


The  first  section,  A,  shows  the  results  of  the 
linear  theory  (Anderson  and  Archambeau,  1964)  using 


J 
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the  exact  velocity  structure.  The  inversion  results 
in  a  highly  oscillating  Qs  structure  even  with  error 
free  data  and  an  exact  velocity  structure.  The 
oscillations  represent  up  to  a  30%  error,  yet  the 
forward  problem  agrees  with  the  exact  QL  in  Figure  1 
to  well  within  0.5%. 

In  section  B  is  the  result  of  simultaneously 

inverting  for  V  and  Q  using  the  exact  treatment, 
s  s 

The  Q  structure  was  perturbed  in  an  oscillating 
s 

manner  towards  the  linear  solution  to  determine 
whether  the  iterations  would  move  in  this  direction 
or  towards  the  actual  model.  It  is  important  to  stress 
that  the  linear  solution  represents  a  global  minimum 
for  that  method  and  a  local  minimum  for  the  exact 
formulation.  Considering  the  results  in  section  B, 
it  appears  that  the  overall  oscillations  in  have 
decreased  and  the  velocity  structure  has  been  approx¬ 
imately  recovered.  This  is  really  quite  satisfactory 


because:  1)  the  inversion  is  extracting  twice  as  many 

parameters  as  for  either  Q  or  V  alone  and,  2)  from 


the  derivative 


(Figure  4)  it  is  apparent  that 


phase  and  amplitude  data  are  certainly  not  independent. 

The  stability  here  is  due  to  the  derivative  since 

any  attempt  to  invert  with  simply  the  derivatives 

diverges  wildly  in  both  V  and  Q0. 

d  vs  d  so 

Section  C  shows  a  standard  V^-only  inversion 
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using  simply  the  analytical  derivative  (Figure 

v  '  s 

2).  Since  its  final  velocity  model  (smallest  SSRC 
for  a  suite  of  iterations)  is  not  as  close  as  the  full 
inversion  of  section  B  for  the  same  starting  velocity 
model,  the  amplitude  data  actually  adds  resolution  and 
stability  to  the  velocity  inversion. 

An  obvious  result  that  emerges  from  this  study  is 
that  attempting  to  resolve  a  number  of  layers  with  data 
of  limited  bandwidth  can  result  in  an  ill-conditioned 
system.  This  fact,  coupled  with  an  approximate  kernel, 
leads  to  an  oscillating  solution  which  is  far  from  the 
exact  model.  The  exact  formulation  leads  to  a  result 
which  is  closer  to  the  actual  model,  but  the  con¬ 
vergence  is  slow.  One  approach  to  this  problem  is 
to  use  more  sophisticated  inversion  schemes  (Der 
et  al . ,  1970;  Jackson,  1972;  Wiggins,  1972)  which 
generally  achieve  better  stability  at  the  price  of 
decreased  resolution. 


CONCLUSION 

An  exact  formulation  is  used  in  calculating  the 
dispersion  and  attenuation  for  Love  waves  in  a  layered 
linear  viscoelastic  half-space.  The  exact  method  is 
used  to  demonstrate  the  instability  inherent  in  using 
the  approximate  linear  inversion  kernel  (Anderson  and 
Archambeau,  1964).  Although  the  approximate  theory 
agrees  well  within  0.5 %  of  the  exact,  in  a  forward 


sense,  the  inversion  can  result  in  alternating  high  and 
low  Q  layers  where  none  existed.  The  reason  for  this 
instability  is  found  in  the  neglect  of  the  derivative 
.  It  is  thought  that  this  result  strongly  favors 
simultaneous  inversions  for  both  velocity  and  atten¬ 
uation  of  surface  waves  which  may  result  in  better 
resolution  for  both  velocity  and  attenuation.  The  use 
of  dispersion  data  only  may  be  the  source  of  discrep¬ 
ancies  in  the  inversion  of  surface  wave  data  (McEvilly, 
1964)  which  are  sometimes  interpreted  in  terms  of 
anisotropy.  Studies  currently  in  progress  involve  the 
simultaneous  inversion  of  both  Love  and  Rayleigh  wave 
data  using  exact  kernels  and  hopefully  will  further 
resolve  this  matter. 

In  applying  this  formulation  to  observational  data 
two  points  must  be  emphasized.  First  the  phase  data 
must  be  corrected  for  anelastic  dispersion  as  dis¬ 
cussed  by  Kanamori  and  Anderson  (1977)*  One  approach 
is  to  correct  all  of  the  data  to  a  convenient  reference 
frequency  before  applying  the  inverse  method.  The 
second  consideration  is  the  effects  on  the  velocity 
and  structure  due  to  errors  in  both  phase  and 

0 

amplitude  measurements.  The  most  general  conclusion 
is  that  the  resolution  will  be  reduced  (Wiggins,  1972). 
However,  since  the  relative  errors  in  phase  and  amp¬ 
litude  measurements  are  very  different,  the  coupling 
between  the  errors  and  the  model  parameter  requires  a 


/ 


11? 

very  detailed  treatment  (Der,  1972),  and  is  not  within 
the  scope  of  the  present  paper.  In  general  though, 
the  magnitude  of  the  ^  y  term  mandates  accurate 
amplitude  data  if  it  is  to  contribute  significant 
information  to  the  inversion  process. 
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TABLE  1 


A  standard  model  for  a  crust  over  a  mantle 


(Harkrider,  1968). 

Identical 

layers 

were  introduced 

into  the  crust  to 

permit  greater  freedom  in  the 

inversion  process. 

Vg(km/sec ) 

(cgs) 

Qs 

thickness (km) 

3.60 

2.80 

100 

10 

3.60 

2.80 

100 

10 

3.60 

2.80 

100 

10 

3.60 

2.80 

100 

10 

4.50 

3.30 

100 

Oo 

v^r\  rvi  *  w\  ro 


TABLE  2 


Results  of  inversion  methods  using  data  generated 
from  model  of  'fable  Is  A)  linear  theory  inversion  of 


Q.  onlys  B)  exact  theory  inversion  for  <2  and  V  j  C) 

5)  S3 


inversion  for  V  using  only  the  derivative 


H  *3 


(SSRQl, 


sum  square  residuals  Q^s  SSRC,  sum  square  residuals 
c;  NI,  number  of  iterations.) 


Layer 

Qs 

Initial 

V 

3 

2SRQ 

SSRC 

^s_ 

Vs 

Final 

SSRQ 

SSRC 

A 

X 

100 

3.60 

122 

3.60 

3.xlO'3 

2.X10"12 

2 

100 

3.60 

76 

3. 60 

3 

100 

3-60 

130 

3. 60 

u 

100 

3.60 

94 

3.60 

5 

100 

6.50 

100 

4.50 

B  1 

1  10 

3-20 

.08 

0.1 

100 

3-59 

9-xlO'10 

l.xlO"9 

2 

90 

3-20 

98 

3-63 

3 

110 

3.20 

101 

3-57 

a 

90 

3.20 

100 

3.62 

5 

100 

6.50 

100 

4.50 

c  l 

100 

3.20 

.0? 

0.1 

100 

3-53 

3-xlO'3 

-4 

l.xlO 

2 

100 

3.20 

100 

3.46 

3 

100 

3.20 

100 

4.28 

4 

100 

3.20 

100 

2.98 

5 

100 

6.50 

100 

4.50 

ive 
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II.  INVERSION  OF  RAYLEIGH  WAVE  DATA  FOR  VELOCITY  AND 


ANELASTICITY  USING  EXACT  KERNELS 

The  following  treatment  will  be  precisely  that  of 

the  previous  except  Rayleigh  waves  will  be  considered. 

The  formulation  again  is  that  of  the  Haskell-Thompson 

method  extended  to  include  anelastic  attenuation 

(Silva,  1976b)  and  applied  to  Rayleigh  waves.  The  phase 

quality  factor  (Qr)  is  defined  exactly  as  in  section  I. 

An  additional  quality  factor  is  introduced  here  and 

is  termed  the  group  quality  factor.  It  is  defined  by 

the  following  equation  C  QR  =  U  (Brune,  1962).  It 

was  neglected  in  the  Love  wave  treatment  because  for 

a  constant  Qg  structure  ^  for  Love  waves  is  frequency 

independent  and  Qy  =  Qg.  This  holds  for  Rayleigh  waves 

only  if  Q  =  Q  =  constant  and  then  Q,  =  (2  =  Q  ,  The 

s  p  ups 

structure  considered  is  identical  to  that  of  the  pre¬ 
vious  treatment  (Table  1)  with  a  constant  Q  =  200 
structure  added  for  P-wave  attenuation. 


In  Figure  1  is  shown  C,  U,  and  QR  for  the 
crustal  model  and  represents  the  error  free  data.  The 
group  velocity,  U,  is  calculated  using  a  variational 
method  applied  to  dissipative  media  (Appendix  3). 


In  Figure  2  is  shown  the  modulus  and  phas^  of  the 
surface  displacement  ratio  ^0  for  the  attenuating 

W0 

model.  The  phase  is  not  very  different  from  the 
elastic  (90°)  and  the  elastic  modulus  is  within  1 % 
of  the  attenuating  model. 
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INVERSION 


As  in  the  Love  wave  inversion,  the  Taylor  series 
is  truncated  past  the  linear  term  with  C  and  data 
inverted  for  V  and  Q  .  The  phase  quality  factor  is 
favored  as  an  inversion  parameter  rather  than  since 
its  use  is  not  contaminated  by  errors  in  determining 
the  group  velocity.  The  derivatives  are  calculated 
numerically  using  a  parameter  change.  A  10$ 
parameter  change  was  found  to  be  optimum  in  terms  of 
stability  and  speed  of  convergence;  however,  mode 
tracking  proved  to  be  too  difficult  for  this  structure 
and  period  range  and  therefore  mandated  a  smaller 
perturbation  (plots  are  for  a  10$  change). 

Considering  now  the  plots  of  the  derivatives 
(Figures  3-1*0  we  can  note  some  interesting  features 
which  give  insight  into  the  controlling  features  of 
surface  wave  inversions.  In  Figure  3  is  shown  the 

£?C 

derivative  ?  .  The  significant  features  are  the 

/V  s 

degree  of  independence  of  the  curves  and  the  presence 
of  peaks.  Comparing  these  derivatives  with  those  for 
Love  waves  (section  A,  Figure  2)  it  is  apparent  that 
the  Rayleigh  phase  velocity  has  much  more  specific 
information  with  respect  to  material  shear  velocity 
for  this  structure  and  period  range.  This  observation 
is  substantiated  by  comparing  the  shear  velocity  only 
inversions  for  Love  and  for  Rayleigh  waves  (Love; 
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section  1,  Table  2;  Rayleigh;  Table  1). 


Layer  Love 

1  3*53  km/sec 

2  3.46 

3  4.28 

4  2.98 

5  4.50 


Exact  Model 

3.60  km/sec 

3.60 

3.60 

3.60 
4.50 


Rayleigh* 

3.60  km/sec 

3.60 

3.60 

3.60 
4.50 


It  is  quite  obvious  then  that  in  this  case  the  shape 
of  the  derivatives  give  an  excellent  indication  of 
inversion  effectiveness.  A  look  at  the  next  plot 


( 


,  Figure  4)  indicates  that  considerable  dif- 


,)c 

ficulty  would  be  encountered  in  trying  to  extract  the 
four  layer  P-wave  velocities.  Here  a  reduction  in 
layering  is  in  order  and  one  can  be  guided  in  layer 
reduction  by  the  curves.  However,  the  magnitude  is 


also  of  importance  because  the 


jL 


competing  with  the 


2 


<)Vt 


derivatives  are 


derivatives  and  are  roughly 


an  order  of  magnitude  smaller.  This  would  result  in 
putting  most  of  the  adjustment  in  the  shear  velocity 
structure  and  therefore  excludes  simultaneous  deter¬ 
mination  of  both  P  and  SV  velocity  structures. 

In  Figure  5  is  shown  the  first  of  the  attenuation 
derivatives.  As  expected,  it  is  small  (0(10"^)) 
indicating  a  second  order  effect  of  attenuat  on 
phase  velocity.  An  interesting  feature  is  the 


*  Extracted  from  given  tables. 
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appearance  now  of  some  double  peaks  and  greatly  in¬ 
creased  structure.  However,  the  peaks  here  are  less 

separated  than  those  of  ^  £  (Figure  3)*  Also, 
jj  t)r  s  p  q 

comparing  —  (Figure  5)  with  (Figure  9)  it 

is  apparent  that  the  relationship  (Appendix  4) 


<K  „  .  4  °R2  9s2  Jo 

Wl  c  Vs  ^  Qs 


holds. 

The  next  plot,  (Figure  6),  again  shows  the 

lack  of  specificity  of  Rayleigh  waves  to  P-wave 
parameters.  We  may  again  invoke  the  previous  relation¬ 
ship  and  apply  it  to  the  P-wave  parameters  citing 

)  Q 

r:,  (Figure  10) : 

r  p 


£  qr 

T^Tp 


4  Q  ^  Q  ^ 
^  WR  wp 


In  Figure  7, 


P 

(?  qR 


it: 


,  we  begin  with  the  Qt 


atives.  The  similarity  of  this  plot  with 
3)  is  predicted  by  the  relation  (Appendix  4) 


deriv- 

(Figure 


and  by  an  approximate  theory  of  surface  wave  prop¬ 
agation  in  anelastic  media  (Anderson  and  Archambeau, 
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curves  indicates  that  with  the 


1964;  where  now  C  and  V  are  elastic  velocities  and 

s 

results  in  a  linear  inversion  for  attenuation  parameters). 

Jq  r 

The  shape  of  these  \  q1 

correct  velocity  structure  and  correct  Q  structure, 
reasonable  success  should  be  expected  in  a  Qg  inver¬ 
sion.  The  correct  V^,  Vg,  and  Qp  structure  is  required 
in  order  that  the  correct  QR  is  obtained  in  derivative 


calculations . 

The  next  plot,  ^ Q*  (Figure  8),  again  shows 
relatively  little  structure  and  is  indicative  of  P 
parameter  insensitivity  in  Rayleigh  waves.  It  can  be 


QR 

P 


compared  with 


ITP 


(Figure  4) . 


t)  Qr 


Figure  9  shows  the  derivative  fys  •  Again,  the 
double  peak  structure  is  present  (as  earlier  compared 
with  q  in  Figure  5)  and  we  have  increased  structure 

at  the  cost  of  decreased  independence.  The  final  QR 

D  qr  k 

derivative,  p  (Figure  10),  was  earlier  compared 

with  Jq  (Figure  6)  and  further  shows  little  hope  of  a 

*  P  J  Qr 

V  inversion  simultaneously  with  V  ,  as  — -  is  an 


order  of  magnitude  smaller  than 


.  s' 

t)  Qr 

v  V"s 


and  shows  much 


less  independence .  This  comparison  is  analogous  to 


the  iih  and 


comparison.  An  interesting  idea 


might  be  a  mixed  inversion.  That  is,  use  C  for  V 


and  Qr  for  Vp  since 


P  Qr 


appears  to  be  less  linearly 


Jc  "  v  <>Qr 

dependent  than  ^  and  conversely  appears  to 

be  more  linearly  dependent  than  •  Future  work 

along  these  lines  is  needed. 


The  final  plots  (Figures  11-14)  show  the  change 

in  modulus  and  phase  of  the  surface  displacement  ratios 

Uq/Wq  for  a  change  in  Vg  and  Qg.  The  derivatives  show 

that  both  the  modulus  and  phase  are  much  more  sensitive 

to  the  velocity  structure  than  the  Q  structure  for 

s 

this  model.  They  also  indicate  that  neither  con¬ 
tributes  information  independent  of  C  or  Q^. 

TESTS  WITH  SYNTHETIC  DATA 

As  in  the  previous  treatment  of  Love  wu^es,  the 
40  km  thick  continental  crust  was  subdivided  into 
four  identical  layers  (section  1,  Table  1).  Both 
the  crust  and  uppermantle  were  assumed  Poisson  solids 
for  P-wave  velocities  which  are  kept  fixed.  Figure 
1  shows  the  synthetic  data  (C  and  QR)  with  the  same 
period  range  and  number  of  data  as  in  the  Love  wave 
inversion  (120  to  5  seconds  with  5  second  intervals). 

In  Table  1  are  shown  the  results  of  both  the 
linearized  inversion  (Anderson  and  Archambeau,  1964) 
and  the  exact  theory  (both  use  singular  value  de¬ 
composition  with  the  non-linear  appropriately  scaled 
(Wiggins,  1972)). 

The  first  section,  A,  shows  the  results  of  the 
linear  inversion  for  Q  using  the  exact  velocity 
structure.  The  Q  structure  is  far  from  the  true 

o 

structure  (Q„  =  100  sill  layers)  and  the  oscillations 
s 


/ 


are  similar  to  linearized  Love  wave  inversion  (Table  2, 

section  A).  Since  the  error  here,  as  much  as  39$,  is 

greater  than  the  linearized  Love  wave  inversion,  part 

of  the  blame  must  be  put  on  the  approximate  data 

kernels.  That  is,  the  partial  derivatives  used  in  the 

linearized  Q  inversion  are  numerically  calculated  for 
s 

Rayleigh  waves  and  are  analytical  derivatives  for  Love 
waves.  The  numerical  derivatives  in  the  linear 
Rayleigh  formulation  yield  approximately  a  5$  error 
in  the  forward  calculation.  However,  the  error  may 
be  significantly  larger  for  inverse  operations.  A 
simple  calculation  may  serve  to  point  this  out  as  it 
indicates  that  small  errors  in  inversions  kernels  can 
be  magnified  especially  in  large  condition  number 
systems. 

Let  j 

D  =  (L  +  e)M  =  L( I  +  L"1e)M 
where  D  is  some  data  vector,  L  some  data  kernel  with 
error  terms  e,  and  M  the  model.  Now  if  e  is  say  5$ 
of  L  in  some  norm,  then  D  will  have  an  error  of  the 
same  magnitude.  On  the  other  hand,  inverting  the 
system  we  have 

M  =  (I  +  L-1e ) _1  L-1  D. 

Now  if  we  let  M'  =  L-'*'  D  be  the  real  model  (error 
free)  then 

M  =  (I  +  L“1e)"1  M' 
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or 

MJtf  (I  -  L~^e  )M*  ,  |IL“1e  1 1  <1 
then  the  model  we  get  (M)  can  differ  from  the  true 
model  (M*)  by  more  than  5%  even  though  L~^  is  well 
behaved  (i.e.f  ||L~^e  IK1)  •  This  should  discourage 
non-iterative  inversions  with  approximate  kernels. 

This  error  magnification  is  probably  contributing 
to  the  very  poor  inversion  results  of  the  linear 
theory  and  should  further  underscore  the  need  for 
extreme  caution  in  evaluating  its  results.  However, 
the  comparison  between  the  linear  approximation  and 
the  exact  theory  is  not  unjust  because  the  error 
certainly  cannot  account  for  all  of  the  discrepancy 
and  the  same  numerical  derivatives  are  used  in  the 
exact  formulation  (future  work  will  investigate  the 
effect  of  analytical  derivatives). 

In  section  B  are  shown  the  results  of  the  exact 

theory  inversion  for  V  and  Q_,  with  an  exact  Q_ 

s  s  p 

structure  (held  constant).  The  results  are  quite 
satisfactory  as  both  the  V  and  Q  structures  have 
been  recovered.  Indeed,  the  results  are  superior 
to  the  Love  wave  inversion  as  was  expected  from  the 
analysis  of  the  derivative  plots.  The  velocity  only 
inversion  is  shown  in  section  C,  shows  good  convergence, 
and  is  far  superior  to  the  Love  wave  velocity  inversion. 
It  appears  then,  that  for  this  model  the  amplitude  data 
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is  not  essential  in  the  velocity  inversion.  Perhaps 
a  further  perturbed  starting  structure  would  require 
the  amplitude  data  to  stabilize  the  velocity  inversion. 
This  was  not  attempted  since  more  than  a  10  percent 
change  is  pushing  the  Taylor  approximation  and  it 
was  desirable  to  keep  the  same  starting  models  for 
both  Rayleigh  and  Love  wave  inversion. 

In  order  to  investigate  the  effect  of  Qp  on  the 
inversion,  it  was  fixed  at  the  wrong  value  (250)  and 
the  V  and  Q  inversion  repeated.  The  results  are 

S  S 

shown  in  section  D  and  are  somewhat  surprising.  One 

might  naively  assume  the  effect  of  Qp  to  be  small 

since  it  is  twice  Q  .  However,  as  section  D  shows, 

s 

fixing  it  at  the  wrong  value  (and  with  only  a  25 

percent  discrepancy)  has  a  disasterous  effect  on  Q 

i )  Qr 

and  affects  Vg  through  •  In  light  of  these 

results,  Qp  was  allowed  to  float  using  the  same 

starting  model  as  section  D.  The  hope  is  that  this 

will  release  the  inversion  and  allow  it  to  converge 

to  a  closer  Q  and  V  structure.  The  results  are 
s  s 

shown  in  section  E  and  indicate  that  this  is  indeed 


the  case  as  both  the  Qg  and  Vs  structures  have  essent¬ 
ially  been  recovered  while  the  crustal  Qp  has  moved 
the  correct  way.  An  attempt  at  using  the  linear 


theory  for  a  simultaneous  inversion  for  Q0  and  Qp 
yielded  extremely  poor  values  for  both  and  is  not 


shown. 
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CONCLUSION 

Tests  with  synthetic  data  in  inverting  Rayleigh 
wave  phase  and  attenuation  data  for  velocity  and 
anelasticity  indicate  a  greater  stability  and  reli¬ 
ability  for  an  exact  formulation  over  the  appropriate 
linear  formulation  (Anderson  and  Archambeau,  1964). 

In  addition,  when  inverting  for  anelasticity  Qp  has 
a  first  order  effect  and  must  be  considered  as  an 

inversion  parameter  along  with  Q  . 

s 


f 
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TABLE  1 


Results  of  inversion  methods  using  data  generated 

from  model  of  Table  1  (Section  A)i  A)  linear  theory 

inversion  for  Qg  only;  B)  exact  theory  inversion  for 

Q  and  V  ;  C)  inversion  for  V  using  only  the  deriv- 
s  s  s 

ative  ;  D)  inversion  for  V_  and  Q„  with  wrong  Q_i 
4  vs  s  s  p 

E)  inversion  for  Vgl  Qg,  and  Qp.  (SSRQ,  sum  square 
residuals  Q^:  SSRC,  sum  square  residuals  Cs  Nl, 
number  of  iterations  in  a  suite  of  20.) 


Initial 


Final 


Layer  Qp 

Vs 

SSRQ 

SSRC 

Vs 

SSRQ 

SSRC 

NI 

A 

1 

3.60 

88 

3.60 

2 

3-60 

S3 

3.60 

3 

3.60 

120 

3.60 

0 

4 

3.60 

61 

3.60 

5 

4.50 

86 

4.50 

B 

1 

200 

no 

3.20 

200 

100 

3.6O 

2 

3 

200 

200 

90 

no 

3.20 

3.20 

l.xlO'7 

?.xl0'2 

200 

200 

99 

100 

3.60 

3.60 

3-xlO 

-*2.xl0-10 

7 

4 

2C0 

90 

3.20 

200 

99 

3.60 

5 

210 

100 

4.50 

200 

100 

4.50 

C 

1 

3.20 

3.60 

2 

3 

3.20 

3.20 

7-xlO*2 

3.60 

3.60 

2.xl0"10 

6 

4 

3.20 

3.60 

5 

4.50 

4.50 

D 

1 

250  no  3.20 

2  50 

98 

3.61 

2 

250  90  3.29  ,  , 

250  no  3.20  a.xio"-  7.xio'<: 

250 

98 

3-55  _6 

-6 

3 

250 

104 

3.69  9-xlO 

3-xlO  0 

6 

4 

250  90  3.20 

250 

90 

3.54 

5 

250  100  4.50 

250 

98 

4.50 

E 

1 

223 

99 

3.60 

2 

227 

99 

3-59  a 

_  q 

3 

Same  as  (D) 

229 

99 

3.60  l.xlO 

7.xl0  y 

13 

4 

22-- 

97 

3.59 

5 

274 

98 

4.50 

4*4 


indicated  with  layer  5  "the  half-space. 
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Summary 

A  method  is  presented  which  allows  the  variational  formulation  for 
elastic  surface  waves  to  be  extended  to  the  case  of  dissipative  media.  With 
this  formulation,  correct  to  second  order  in  the  loss,  Rayleigh’s  principle 
can  be  applied  to  perturbations  of  the  Rayleigh  quotient  to  yield  group 
velocity  without  numerical  differentiation.  Other  perturbations  can  be  used 
to  find  the  change  in  phase  or  group  velocity  due  to  changes  in  loss, 
density,  or  moduli. 


Introduction 

The  use  of  Rayleigh's  principle  is  of  considerable  importance  in  surface-wave 
calculations.  As  suggested  by  Meissner  (1926)  and  Jeffreys  (1959, 1961)  and  amended 
by  Harkrider  (1968)  it  has  replaced  numerical  methods  with  an  exact  formulation  for 
calculating  group  velocities.  In  addition,  Rayleigh’s  principle  may  be  used  to 
calculate  the  effect  on  the  velocity  dispersion  due  to  small  perturbations  in  the  elastic 
parameters.  This  type  of  information  is,  of  course,  most  essential  in  solving  inverse 
problems. 

However,  in  its  usual  form,  Rayleigh’s  principle  is  inadequate  for  dissipative 
systems  since  its  use  requires  equating  the  time  average  kinetic  and  potential  energies. 
This  result  follows  from  the  virial  theorem  and  imposes  a  vanishing  Lagrangian  for 
the  system. 

The  purpose  of  the  following  development  is  to  construct  a  Lagrangian  for  a 
non-conservative  system  which  vanishes  and  then  apply  Rayleigh’s  principle  to  the 
resulting  equation.  This  can  be  done  by  writing  a  Lagrangian  for  two  systems 
operating  simultaneously  with  one  losing  energy  as  the  other  gains  energy,  so  that  the 
total  energy  is  conserved.  An  estimate  of  the  error  introduced  in  combining  the 
systems  is  given  in  the  Appendix. 

The  development  will  be  for  Love  waves  in  order  to  simplify  the  equations.  The 
extension  to  Rayleigh  waves  follows  naturally. 


Formulatioi] 

As  a  working  example,  we  shall  consider  Love  waves  propagated  along  the 

445 


t 


I 


A 
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surface  of  a  vertically  heterogenous  (layered)  attenuating  half  space.  The  usual 
equation  to  be  satisfied  in  each  layer  is 


V2u+  —  V2ii  =  pit  (1) 

where  u  is  the  transverse  particle  displacement,  p  is  the  medium  density,  and  /i„  and 
//,  are  the  conservative  and  non-conservative  Lame  parameters,  respectively,  which 
are  in  general  frequency  dependent  (Borcherdt  1971;  Silva  1976).  We  take  as  the 
solution  to  equation  (I) 


u  =  v(z)  exp  ( -  Ax  x)  cos  (Pxx—(ot)  (2) 

where  Px  is  the  spatial  frequency  such  that  the  horizontal  anelastic  phase  velocity  c 
is  given  by  c  —  a>l  Px  and  Ax  is  the  horizontal  attenuation  factor.  Application  of  the 
usual  boundary  conditions  results  in  the  continuity  of  both  Px  and  Ax.  They  then 
represent  eigenvalues  to  be  determined  for  propagating  modes.  The  attenuation  factor 
may  be  defined  as 


_  iH  r~1+v^t+Qi-~a)l*  02 

*  c  I+VO+Gl"2)  *  2 cQl 


(3) 


where  QL  represents  the  effective  quality  factor  for  the  spatial  decay  of  the  surface 
wave.  Equation  (3)  can  then  be  interpreted  as  the  projection  of  the  attenuation  in  the 
propagation  direction.  The  shear  quality  factor  Qs  is  defined  as 


In  E  pR 


(4) 


where  A  £  is  the  energy  lost  and  E  the  peak  energy  stored,  both  per  cycle. 

At  this  point  we  introduce  the  mirror  image  system  which  is  taken  to  exist 
simultaneously  with  the  original  field.  Since  this  field  gains  energy  exactly  as  the  first 
loses  energy,  we  can  write  the  field  equation  as 


/<»  W-  —  V”«  +  =  pit* 

O) 


(5) 


with  the  corresponding  solution, 

u  =  r(z)  exp  (  >,  x)  cos  ( Px  x  -  w/).  (6) 

It  is  seen  then  that  the  only  difference  in  the  systems  is  in  the  sign  of  p,.  This  then 
affects  only  the  sign  of  the  attenuation  constant  and  yields  the  same  spatial 
frequencies  as  equation  (2). 

In  order  to  construct  a  conservative  Lagrangian  density  we  make  use  of  the 
following  form  (Morse  &  Feshbach  1953;  Moisewitsch  1966). 


,  |i|  ru+  till  cu+  I 

X  =  puu+-±—  V*fuu+-uu  ]~pK  U--T-  +  T-r-  • 
ot  [f.x  r.v  cz  J 


(7) 


This  yields  an  invariant  Lagrangian,  since  m+  increases  in  amplitude  exactly  as  u 
decreases. 
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The  corresponding  Euler- Lagrange  equations  become: 

_  d_d&  _  d_d&  _  _  Q 

Du  dt  Du  dx  dux  dz  Du. 


d_  D&_ 
dt  du  + 


d  D<£  d  DUP 


dxdux+  dz  du. + 


Direct  substitution  of  equation  (7)  into  equation  (8)  yields  equations  (1)  and  (5). 
We  now  substitute  the  solutions  (equations  (2)  and  (6))  into  equation  (7)  and  then  space- 
time  integrate  the  Lagrangian  density  to  yield  the  Lagrangian  of  the  combined  system. 
In  the  time  and  x  integrations,  where  periodic  solutions  were  assumed,  the  integrations 
result  simply  in  averages. 

Defining  L  as  the  Lagrangian  of  the  combined  system,  we  can  write: 

L=  to2  j  pvv+  dz—(Px2  —  Ax2)j  ftRvv+  dz—  J  /iR  v'v+'dz .  (9) 

The  first  term  is  interpreted  as  a  time  average  kinetic  energy  and  the  remaining 
two  terms  as  a  time  average  potential  energy.  This  allows  the  application  of  the 
virial  theorem  (Moisewitsch  1966)  which  states  that  for  a  conservative  system  which  is 
quadratic  in  its  potential  energy  the  time  average  potential  and  kinetic  energies  are 
equal. 

We  may  then  put  equation  (9)  onto  the  form  of  the  Rayleigh  quotient: 

<o2I0  =  (Px2-Ax2)I,+I2  (10) 

with  the  energy  integrals 

l0  =  J P«’+  dr;  /,  =  J  //R  w+  dz\  I2  =  J  /<R  v'v+'dz. 


Application  of  Rayleigh's  principle 

Rayleigh’s  principle  states  that  for  a  given  Rayleigh  quotient,  as  in  equation  (10), 
any  eigenvectors,  correct  to  first  order,  will  yield  the  eigenvalue  appropriate  to  that 
mode  correct  to  second  order.  This  can  be  stated  more  precisely  by  considering 
the  perturbation  to  equation  (10)  due  to  a  smalt  change  in  v  and  i,+. 

w2dl0  =  (Px2-Ax2)dIl+5I2  (11) 

where  the  perturbed  energy  integrals  are  given  by: 

<5/„  =  J pd(cv+)dz;  <5/,  =  J ptt5(vv+)dz;  5I2  =  J pR5(v' v+’)dz . 

Making  use  of  equation  (11)  the  new  eigenvalues  corrrect  to  second  order  may  be 
calculated. 

(<o+<MJ[/0+<5/o]  =  [(P.+SP^-iA.+dA^.lh+SI^+lh+SI,).  (12) 

Using  equations  (10)  and  (11)  and  neglecting  second  order  in  small  quantities 

equation  (12)  becomes: 

cod(ol0  =  ( PxSPx-AxSAx)It .  (13) 

and  using  c  =  co/Px  for  phase  velocity,  V  -  5oilSPx  for  group  velocity,  we  can  rewrite 
equation  (13)  in  the  following  form: 


A*6A,) 

PJPA 


(14) 


At  this  point  Rayleigh's  principle  may  again  be  invoked  and  the  elastic  eigenvectors 
may  be  used  in  the  calculation  of  the  energy  integrals. 
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in  order  to  estimate  the  perturbation  in  AX(&AX)  due  to  a  frequency  perturbation,  we 
use  equation  (3)  and  neglect  the  term  containing  8QJ5w.  This  is  correct  to  at  least 
second  order  in  QL~ '  since  QL  is  a  slowly  varying  function  of  frequency  and  the  term  has 
a  coefficient  of  Qi~2.  (It  can,  in  fact,  be  estimated  using  a  difference  scheme  since 
QL(u>)  is  known  from  the  eigenvalues.) 

With  this  approximation  equation  (14)  becomes 


which  reduces  to  the  usual  elastic  expression  as  QL  -+  oo  and  v*  -*  v.  This  expression, 
which  is  an  approximation  correct  to  at  least  second  order  in  the  loss,  is  thought  to  be 
preferable  to  numerical  differentiation  of  the  phase  dispersion  curve. 

f  n  order  to  calculate  the  effect  on  QL  due  to  a  small  change  in  (?s  we  must  first  find  an 
expression  for  /iK  in  terms  of  t\  and  Qs.  This  can  be  done  by  assuming  a  solution  of 
equation  (1)  of  the  form 

u  =  u0  exp(  —  »'(K.X  — «/))  (16) 

where  K  =  P  — /A  (Borcherdt  19/ j).  Substituting  this  into  equation  (1),  and  defining 
the  shear  velocity  v,  as  that  of  homogeneous  waves  (P  parallel  to  A)  we  can  write: 


flR 


v*1  ri+v'(l+es~2)  I 

p  2  1  i+or2  J 


(17) 


By  substituting  equation  (17)  into  I,  and  I2  we  are  able  to  explicitly  calculate  the 
change  in  the  Love  wave  quality  factor  8Q,  due  to  a  change  5QS  in  the  energy  integrals. 
We  can  write  the  entire  perturbed  Rayleigh  quotient  at  constant  frequency  as: 

<’>2[lc+M0]  =  {(Pa+.SP()2-(4(+<\4,)2].(/,+^i+^s/.] 

+  [/j  +<5/j  +<5/js  1 2}  (18) 


where  we  have  defined 


,  .  r  rsJ  * ,  ri+vn+Qs"2)i . 

■v'.  -  ;/’—*<■  <l— 

,  ,  _  f  r,1  ,  ,,,ri+y/(l-fgs~2)l  , 

,Q*  1  -  b-j1  »[  ,+0s-2  Jrf‘- 


(19) 


Using  equations  (10)  and  (11)  and  neglecting  second-order  terms  in  small 
quantities  we  can  reduce  equation  (18)  to: 


0  =  2[PX8PX-AX8AX]I,  +  [PX2~AX2]8Q%1,  +<5qs  12. 
Writing  equation  (3)  as 

-1+v/(1+0l'2)1’  j 


a,*  PJlQiYJ(Qi) 


(20) 

(21) 


l+s/d  +Ql 

we  can  express  the  perturbation  in  Ax  due  to  a  change  in  gL: 

t>Ax=bPxf(Q0  +  Px8f(Q0.  (22) 

Neglecting  the  first  term  since  it  is  at  least  second  order  in  the  loss  we  then  have  from 
equation  (20); 

a  a  1 

•  -y  =  HP'2-AS)S*l,  +SQshl 


(23) 
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and  to  the  same  order  we  have  from  equation  (3) 

6Qux  -2QL26f(QL).  (24) 

This  formulation  then  allows  estimates  to  be  made  on  the  change  in  the  Love-wave 
quality  factor  or  attenuation  factor  due  to  a  small  change  in  the  shear  quality  factor. 
Further  perturbations  involving  rs  and  p  can  be  made  along  similar  lines  (Harkrider 
1968). 

The  extension  to  Rayleigh  waves  will,  of  course,  be  more  tedious  since  equation  (9) 
w  ill  now  contain  four  integrals.  In  addition,  equation  (18)  now  must  consider  Qp  as 
well  as  Qs.  The  situation  is  unpleasant  though  tractable  and  the  group  velocity  is 
presently  being  programmed  by  the  author. 
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Appendix 


In  order  to  estimate  the  error  introduced  by  the  formulation  presented  in  this  paper 
we  shall  consider  a  plane  shear  wave  propagating  along  the  x  direction  in  a  Voigt 
solid  (Kolsky  1963). 

The  equation  of  motion  can  be  written 


with  the  solution. 


02  if  , 


Px1 


=  pw 


tv  =  exp  ( -  Ax)  cos  (lot—  Px). 


(1) 

(2) 
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We  may  define  the  quality  factor  for  this  system  as: 


Q"'  = 


P 


(M 


Substituting  equations  (3)  and  (2)  into  (1)  wc  arrive  at: 
P2-A 2  = 


a  2p 


~2\  ‘ 


(4) 


On  the  other  hand,  if  we  had  applied  the  Lagrangian  formulation  presented  in  this 
paper  the  results  would  have  been: 


p 


(5) 


which  is  correct  to  second  order  in  the  loss. 
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APPENDIX  3 

DEVELOPMENT  OF  ENERGY  INTEGRALS  AND  GROUP 
VELOCITY  FORMULATION  FOR  RAYLEIGH  WAVES 

Applying  the  formulation  presented  in  Appendix  2 
to  Rayleigh  waves  and  using  elastic  eigenvectors,  the 
energy  equation  becomes 

"*1.  *  (r,‘-  ♦  a  M,  *  xT  ; 

where 

r's  A*  Ifcl‘]‘/Z 

:  *  fe)(  *]  * 

r- ([«■->  fer-  '-tef]-" 

with  U1  and  ^3  the  horizontal  and  vertical  particle 
Wo  Wo 

displacements  normalized  to  the  vertical  surface 
displacement.  The  Lame  parameter  R  is  defined  in 
Appendix  2  with  \  R  similarly  derived  and  fi  is  the 


/ 
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phase  difference  between  the  vertical  and  horii  jntal 
attenuated  displacements. 

Performing  the  appropriate  perturbations  to 
obtain  the  group  velocity  the  Rayleigh  wave  equation 
analogous  to  eq.  (13 )  of  Appendix  2  becomes 

l x.  =  (i  -  fit  V  r. 

Mp»  V  P<IP<J  Px 


The  group  velocity  is  then 


1 


Since  the  elastic  eigenfunctions  are  used  in 
the  energy  integrals  the  integrations  can,  in  principle, 
be  done  analytically.  However,  due  to  the  tedious 
algebra  required  it  was  decided  to  do  the  layer 
integrations  numerically  while  doing  only  the  half¬ 
space  contribution  analytically.  The  integration 
method  adopted  was  that  of  Gauss-Legendre  quadrature. 

The  scheme  uses  the  two  point  computation  applied  to 
successive  points  within  each  layer.  Estimates  of  the 
accuracy  are  routinely  obtained  by  checking  the  balance 
of  the  elastic  energy  equation.  Typical  figures  are 
4-6  function  evaluations  in  a  10-20  km  thick  layer 

_  p 

to  achieve  10  percent  accuracy.  Liquid  layers  are 
included  by  using  the  liquid  layer  matrix  Dm 
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(Chapter  1,  Section  2)  where  appropriate  and  by 
Betting  /:  =  o  in  the  energy  Integrals. 
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APPENDIX  4 

ANALYSIS  OF  INVERSION  DERIVATIVES 


The  purpose  of  the  following  development  is  to 
demonstrate,  analytically,  the  relationships  (mag¬ 
nitude  and  sign)  which  must  exist  between  the  inversion 
derivatives.  The  relationships  follow  from  the  complex 
surface  wave  velocity  being  an  analytic  function  of 
the  layer  complex  velocity. 

Let 


iC0 

C  =  C0  +  20^  • 


iVC 

V  =  V0  +  2Q[ 


(Al) 


where  C  is  the  complex  surface  wave  phase  velocity 
and  V  is  the  complex  shear  wave  velocity  and  both 
have  an  associated  quality  factor.  The  above  expres¬ 
sions  are  actually  low-loss  approximations  and  are  used 
for  ease  of  computation. 


At  this  point  the  Cauchy-Reimann  conditions  are 

»&)  • 

with  respect  to  the  quality  factor  only  (velocity  held 


applied  and  the  partials 
with  respect  to  the  qual 
constant).  The  first  condition  yields; 


are  taken 


/ 
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Applying  the  second  condition  yields* 


resulting  in 

i  aL  4  aL2  9s2  !>  °o 

7^  =  '  co  vo 


(A3) 


These  relations  are  most  useful  in  assessing  the 
relative  magnitudes  of  the  elements  of  the  Jacobian. 
Rewriting  the  Taylor  approximation  (eq.  3)  with 
normalized  parameters  we  have  in  matrix  form* 


Then,  reducing  the  Jacobian  to  two  derivatives  using 
equations  (A2)  and  (A3)  results  in  the  system: 


K 

«L  Vs  JCL  Vs  ‘)aL 

ql 

"s  °L  7^  9L  Va 

Qs 

SCL 

1  Vs  K  vs  ^  CL 

*Vs 

CL 

aL  7^  CL  vsj 

Vs 
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For  the  model  in  Table  1  and  the  corresponding 

partial  derivatives  (Figures  2-5)  the  elements  1,1, 

1,2,  and  2,2  of  the  coefficient  matrix  in  equation  (A5) 

are  all  of  the  same  order  of  magnitude  while  element 

2,1,  due  to  the  Q,  Q  factor,  is  several  orders 

JU  s 

smaller.  These  results  imply  that  element  1,2  is 
significant  compared  to  elements  1,1  and  2,2  and 
cannot  be  neglected. 

Note  that  equations  (A2)  and  (A3)  can  be  used  as 
a  check  on  the  accuracy  of  numerical  derivatives.  It 
can  also  be  easily  shown  that  writing  the  Taylor 
series  approximation  for  Q-^  in  lieu  of  Q  leads  to 
precisely  the  same  matrix  equation  as  for  Q  (eq.  A5) . 

In  other  words,  inverting  for  Q  or  inverse  Q  are 
formally  equivalent. 
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APPENDIX  5 
PHASE  SMOOTHING 

Consider  the  calculation  of  phase  velocity  using 
the  following  formula  (Toksoz  and  Ben-Menahem,  1963)* 

C  =  t ( 2 )  -  t(iftX+  T  +  N) 

Where  Ax  is  the  path  between  stations  and  t(l),  t(2) 
are  times  of  subsequent  Fourier  windows.  This  required 
the  calculation  of  the  phase  difference  (A/)  between 
the  two  spectra.  Since  noise,  including  multiple 
arrivals,  is  generally  present  A fi  may  be  poorly 
behaved.  It  therefore,  becomes  necessary  to  devise 
some  suitable  smoothing  procedure.  The  idea  is  to 
use  some  method  which  is  easily  controllable  so  the 
smoothing  is  not  excessive  and  useful  information 
lost.  The  method  used  here  is  that  of  simply  fitting  A/ 
with  a  sine  series  and  then  retaining  only  those 
coefficients  which  have  the  largest  magnitudes.  The 
coefficients  are  then  modified  by  the  Lancsoz  sigma 
factors  to  reduce  the  associated  Gibbs  phenomena. 
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